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ABSTRACT 

A two-dimensional numerical model is developed to 
predict the temperature and wind fields within a small 
axially symmetric valley cross-section. The thermodynamic 
equation, coupled with the vorticity equation, is solved 
using finite-differences subject to either a surface cooling 
rate which is constant or else dependent on radiative 
transfer processes, following Brunt (1934). The cooling 
rates were obtained from experiments performed in 1977 - 
1978 within the North Saskatchewan River valley in Edmonton, 
Alberta, Ae uae clear sky and light wind inversion 
conditions in which downslope winds developed. Advection, 
turbulent diffusion and radiation processes are incorporated 
into the heat transfer equation ; K-theory 1s used to 
approximate.the latter two. Integrations are performed using 
the Duftort-Frankel formulation for the diffusion terms but, 
due to inherent numerical problems, a forward-time, 
centered-space method is proved Superior in this 
application. A forward-time, upstream-Space formulation is 
used to estimate the advection terms. Liebmann sequential 
relaxation is used to determine the streamfunction field. 

The model predicts attainment of quasi-steady-state 
conditions within about 20 minutes of integration time with 
the development of downslope wind speeds of = 0.8 m s~'. The 
relative contributions which advective, diffusive, and 
radiative processes make are illustrated in different 


conditions and various valley locations. In the early stages 
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(after 200 seconds), radiative transfer dominates at low 
levels. As the slope winds develop, advection processes 
become increasingly important and, after 1600 seconds, cool- 
ing rates due to advection of -4 x 10-* C° s-' are predicted 
4 m above the slope, about twice those due to diffusion or 
radiation. However, when only radiative transfer is 
permitted, which is simulated in this study as a diffusive 
process, it iS able to account for most of the observed tem- 
perature changes alone at low levels. The predicted 
low-level temperature variations with height and the maximum 
Slope winds agree extremely well with the observations. 
Sensitivity of the results to the location of the upper 
boundary of the model and to the presence or absence of a 


flat valley bottom is small close to the slope. 
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Chapter 1 


INTRODUCTION 


1.1 Mountain and Valley Winds 

The phenomenon of mountain and valley winds has been 
investigated by many authors in the past. Wagner, during the 
1930's, published a number of papers giving fairly extensive 
explanations of the theory of such local circulations. Based 
on observations conducted within the valleys of the Alps, he 
concluded that a pressure gradient was created between air 
over the inclined mountain slope and air at the same height 
over the center of the valley as a result of differences in 
the temperature regimes at these two locations. He 
postulated the existence of a double vortex type of circula- 
tion, with weak downslope or drainage winds, ascent above 
the valley floor, and a generally downvalley wind close to 
the valley floor. This latter component was fed partly from 
the side slopes, and also resulted partly from the pressure 
difference between the mountain and the plain. 

Defant (1951) schematically illustrated the successive 
stages in the diurnal cycle of the mountain-valley wind 
system. This is reproduced in Figure 1.1. In (a), at 
Sunrise, upSlope winds begin, but the valley is still colder 
than the plains above, and the downvalley winds are still 
blowing, teduby. the return circulatrzon trom the slopes. it 
Soon Gies Out With further heating. In (pb); the forenoon, 


Slope winds are strong, and a transition from downvalley to 
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upvalley wind is occurring with the valley and plains tem- 
perature approximately equal. In (c), at noon or early after- 
noon, the slope winds are diminishing, the upvalley wind 
is fully developed, and the valley becomes warmer than the 
plains. By (d), late afternoon, the slope winds cease, the 
upvalley wind continues, and the valley remains warmer than 
the plains. During (e), evening, downslope winds begin, the 
upvalley winds diminish, and the valley remains only slightly 
warmer than the plains. By (f), early night, the downslope 
winds become fully developed, with the valley at the same 
temperature as the plains. During (g), the middle of the 
night, the downvalley wind is fully developed, with the valley 
being colder than the teins. ‘the downslope valley wind 
continues. In (h), late night, the downslope wind has 
ceased, with the downvalley wind filling the valley 
completely, and the valley remaining colder than the plains. 
This progression can be interrupted or altered under 
various conditions. Development of cloud cover tends to slow 
the cooling experienced at the ground and greatly weakens 
the valley inversion. This can then result in very low or 
non-existent slope winds. Interaction between the 
large-scale prevailing flow and the thermally induced slope 
wind circulation has been investigated by Tang (1976). He 
showed that a separated cell would form on the lee slope 
during the day, and that one would form on the windward side 
at nighttime. Very strong synoptic scale circulation is able 


to mask the local mountain and valley wind system completely. 
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1.2 Micrometeorology of the N. Saskatchewan River Valley 

The North Saskatchewan River valley bisects the city of 
Edmonton, Alberta (located at 53° 33° latitude, 113° 30° W 
longitude). Typically, it is 50 m in depth and varies from 1 
- 1.5 km across. Along its meandering course (see Figure 
1.2) are both steep-sided valley walls as well as more 
gentle flood plains. Within the city limits, the valley is 
fairly well covered by wooded areas. 

The earliest micrometeorological observations within 
the North Saskatchewan River valley were made in 1958 - 
1959. Klassen (1962) found distinct regimes to exist on the 
plain, in the ravines, and in the main river valley, which 
resembled those suggested by Defant (see Section 1.1). Lower 
surface temperatures and stronger inversions were found to 
occur within the valley at night. Under clear skies and 
light winds (less than 5 m s~'), downslope and downvalley 
winds were generally observed. The speed and depth of the 
wind down the ravines uSually increased to their maximum 
values (1.8 m s~' and the ravine depth, respectively) within 
less than two hours after sundown. Trajectories of fog and 
smoke showed ascending motions above the river and a return 
flow to the valley slopes, in a helical fashion. 

In 1977 and 1978, thirteen field experiments were 
carried See pe Hage (1979) in order to study the influence 
of a small urban valley on air pollutant concentrations. As 
part of the program, air temperatures, winds (and some ver- 


tical profiles of each) were obtained at various valley 
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Figure 1.2 Topography of Edmonton (height contours in 
meters) 
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locations, as well as measurements of carbon monoxide. These 
were supplemented by data from the various rural and urban 
airports in the vicinity. Analyses of the observations, 
especially under conditions of clear skies and light winds, 
clearly suggest the existence of a valley micro-climate 
distinct from that of the city. The general characteristics 
are summarized below. 

In early evening, shortly before sunset, cooling began 
as the radiation received at the ground dropped. Usually, 
the valley became isothermal earlier than the air over the 
city. Inversions developed over both valley slopes and were 
usually more intense than above the urban plains. That slope 
which was shaded first (north- or east-facing slope) 
developed an earlier and stronger inversion (from 8 to 18 
K° (100 m)-'). Differences in temperature of 4 to 6 K° were 
common between the ridge and trough of the valley, with the 
latter location being similar to the rural airport (Hage 
£9992) < 

Downslope winds with typical speeds of 0.4 to 0.6m s7' 
commenced shortly after the inversion formed. Slope wind 
speed was found to be highly correlated with the vertical 
temperature gradient (and corresponding horizontal tempera- 
ture gradient) and must be associated with large pressure 
gradient forces. The slope wind maximum was observed to 
occur about 2 to 4 m above the slope, and, as the inversion 
deepened, it dropped below 0.7 m, the lowest measurement 


level. Above the ridge, a normal increase in wind speed with 
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height was observed. 

A second evening maximum in pollutant (CO) 
concentration was observed over the city and especially 
within. the. valley., This. is.due.to. the development of the 
surface-based inversion. The increased stability produces 
much decreased vertical dilution and can overcompensate for 
the lower pollutant emission rates beyond the rush hour 
traffic peak. Recirculation of air in organized cells may 
also contribute to the high evening concentration within the 


valley. 


1.3 Objectives of the Study 

Almost all previous numerical work in this area has 
been aimed at the description of local mountain-valley wind 
circulations of a physical scale which greatly exceeded the 
dimensions of the valley under consideration. It was 
therfore deemed desirable to investigate the valley circula- 
tion induced by cooling along a relatively shallow slope 
(0.1 - 0.2) of a small urban valley (with a depth of about 
50 m) by the use of a numerical model. The observational 
data already available would be useful not only in model 
verification but also as hints to model development. 

This thesis contributes to a three-part study of the 
micrometeorology of the North Saskatchewan River valley. A 
one-dimensional radiative-conductive model has been 
developed by di Cenzo (1979) to describe the evolution of 


temperature profiles in different parts of the valley. He 
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concluded that the inclusion of advective processes would 
produce a major improvement in the results. This leads 
naturally into the need for prediction of wind fields 
consistent with the surface temperature changes which occur 
as the evening progresses. The thermodynamic variables and 
the corresponding circulation are intimately related and 
only limited success is to be expected if treated 
independently. However, the use of simplified variations in 
the thermal characteristics of the valley can provide some 
insight into these inter-relationships. Rudolph (1980) 
developed a pollutant transport model with an assumed 
double-vortex wind field within a 'V-shaped' valley and 
dispersion of pollutants by turbulent diffusion. Again, in 
order for such a model to predict concentrations accurately, 
the spatial and temporal variations in the thermodynamic and 
wind fields must be incorporated. This thesis describes the 
development of a model capable of predicting such fields 
within a two-dimensional axially symmetric valley cross- 
section. The thermodynamic equation, coupled with the 
vorticity equation, is solved using finite differences 
Subject to either a constant surface cooling rate or one 


which is dependent on radiative transfer processes. 


% 


{ 
| 


=i 


bivow aseeeso7g aviwt 
ebael aid? .esis e1 ada nt 


“a eli 


—, 


Sleit boiw 39° AGi rm be iy 36" 


iwose doidw esposds rueregues « satus , 
al 


18 2esidelyev otmanghanneds auietii serge 


ale: vietsmtont ora acigeiat ih 
costs 
ream” 


ms be? 


betsess $% Bevosqes od ad Sr” 

‘aizayv betiefiqmte to sau erty ‘te 

7 

noe abivorq e469 veliev eng 36 eciveirvetne 
peony 

igqlobus .aqgidenolisiox- Ween? seat 


bemuezs an dtiw febom sxeqensas neat fe 


i Ls 70qgms 
5 2 & 4 £10¢ It 


a F <> 
= onissibexra to sideeas febem a6 


amnye yileixs lenetensmiseowr | 
Pha 
> nbtsevps oimenybonelé ” 
- , 
oeriay Beviog 8: .nez 


[> soetwe Jnevznoo « tenths a 


: 7 ‘ po 
‘stenads3 svitelbhbs? 16 tasbeaage 


Chapter 2 


THE THEORETICAL MODEL 


2.1 Introduction 

This Chapter gives details of the valley wind model 
developed. Assumptions, initial and boundary conditions 
required by the model are outlined. Numerical aspects 
are considered in Chapter 3. 

Ideally, a valley wind model should provide the three 
velocity components and the temperature, pressure and 
density fields naselgty the valley, over it, and above the 
adjacent plains, all as functions of space and time. The 
Initial conditions of geographical location, topography, 
type of ground surface, the radiation away from it as a 
FUNeCETONMOL pOSitwon and time (or insolation ‘onto "the: sur- 
face in the case of daytime heating), and the prevailing 
atmospheric conditions such as the synoptic-scale flow, 
Stability, humidity and cloud presence all influence the 
development of the slope wind in a very complex manner. 
Given this situation, it is desirable to concentrate on the 
most fundamental features of the phenomenon. To this end, 
the model does not possess a detailed boundary layer, nor 
the inclusion of moist processes (other than assumptions 
involving rates of heat loss at the surface given cloud-free 
skies). 

Although computers impose much less of a restriction on 


modelling efforts than at one time, significant increases in 
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running costs and storage amounts would be expected in going 
from two to three space dimensions, or in going from simple 
to real terrain. Checks on the basic method of calculation 
can be made at this simple stage prior to extension into 
more complex regions. As well, the interpretation of results 
generally increases in difficulty as the model complexity 
increases. The likelihood of encountering new numerical 
errors, which are also difficult to identify and interpret, 
increases with complexity. 

With the above simplifications, it is hoped that the 
flow patterns obtained can provide some understanding of the 
transfer processes involved. This understanding can then 
serve aS a basis for analyzing data already available and 
those from future field experiments, as well as for 
understanding more sophisticated models. 

Advection processes are of great importance in the cir- 
culation being studied and, therefore, the complications 
induced by the non-linearity of the equations governing tem- 
perature and velocity changes are necessary. Thus, a numeri- 
cal approach to the problem is required. A finite-difference 
grid model was developed in order to predict the flow and 
thermodynamic fields as a function of space and time. 

The model is initiated with the atmosphere at rest, 
under a variety of stability regimes. A disturbance is 
introduced by cooling the air at the surface and the 


governing equations are solved numerically. 
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2.2 The Governing Equations 

The basic equations used were the equations of motion, 
the thermodynamic equation, and the continuity equation, 
along with the equation of state and Poisson's equation. 


These equations can be expressed, respectively, as 


El Ane Nes ele lat 2 av_ 
ye WORKERS may gk 2G 0+ tis Fe | Fe (2 Deets) 
90 2 e oar 30 
Se DCS Ore [), oa (2.2.2) 
>You 183) ol oo 
=-V°V ae Spr ads (VeV)p 2 Z ca) 
Pe te (ona 
R/C 
z. 2] P (2.2.5) 
PO 


where the last term in (2.2.1) represents the turbulent flux 
of momentum using K-theory in which turbulent transfer of 
momentum is considered in the same manner as molecular dif- 
fusion (which is neglected here). Repeated indices imply a 
summation with respect to that index. Note 6; ;=0 if i#j and 
6;;=1 if i=j. All other symbols have their usual 
meteorological meanings. 

Two-dimensional flow fields were treated in this study 


(in the vertical (z) and cross-valley (y) directions), with 
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the down-valley flow (u) equal to zero. As a result, the 
u-component of (2.2.1) is eliminated. As the velocities to 
be considered are small, the coriolis accelerations were 
negligible compared to other terms in the equations of 
motion. Therefore, axial symmetry was assumed and only the 
fields of one-half the valley cross-section were computed. 
However, differences in topography, slope aspect, etc. must 
be taken into account when comparing the south-facing and 
north-facing slopes and their respective circulation 
regimes. 

As in early treatments of the planetary boundary layer, 
the turbulent flux of momentum was considered to be 
analagous to molecular diffusion, using an eddy viscosity 
coefficient, K,. Just as random motion of molecules over 
some mean free path leads to a transport of momentum on a 
small scale, irregular turbulent motions over some mixing 
length will also cause such a transport. The essential 
Gifference, however, 18 that the molecular coefficients are 
known functions of state whereas the corresponding turbulent 
exchange coefficients depend not only on the characteristics 
of the fluid but especially on the characteristics of the 
flow. Except in a millimeter thin layer at the ground, the 
exchange processes at the molecular level can be effectively 
disregarded. 

In magnitude, the eddy momentum and heat diffusivities 
are much larger than the corresponding molecular values. The 


exact nature of the appropriate eddy values is typically a 
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complicated function of height above the ground, and the 
velocity and temperature profiles. The height above ground 
is important in the definition of an eddy mixing length. 
Close to the ground, the eddy size is physically limited by 
its distance from the ground. Further aloft, a free atmos- 
phere eddy mixing length can be used which is no longer 
constrained by surface friction effects. Mechanical 
turbulence or forced convection is increased when the verti- 
cal wind profile indicates large amounts of shear. As the 
Shear increases, a flow which may have been laminar ini- 
tially will become unstable and develop large eddies which 
can transport heat and momentum quickly. However, when the 
atmosphere is stable, any vertical displacements will be 
Suppressed, and the mechanical turbulence weakened. On the 
other hand, an unstable boundary layer acts to enhance the 
thermal turbulence through buoyancy forces. 

From the turbulent kinetic energy equation, the terms 
describing production of kinetic energy through buoyancy and 
shear can be written as gé’w' /8, and -u’w’ ou/doz, respect- 
ively. All terms have their usual meteorological meanings 
with @) representing the potential temperature of some 
reference atmosphere. The flux Richardson number, Rey is 
defined as the ratio of the buoyancy and shear production 


terms, and iS given by 


t ' 
R, = —&i"—_ (24276) 
os u'w' du/dz 
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if thet heat transfer s upward (erwre>t 0), Riis 


negative because uw < 0 if 00/dz > 0. This implies an 
increase in the turbulent kinetic energy. Upward heat flux 
usually corresponds to 06/az < 0, as in an unstable atmos- 
phere. If the heat transfer is downward (0% win < 0), then Ry 
1S positive and the buoyant production is negative, 
indicating that kinetic energy is lost. Negative values of 
6° w typically correspond to positive values of 06/dz, as in 
a stable atmosphere. If a positive value of R, becomes very 
large, all turbulence becomes suppressed. 

The diurnal changes in solar radiation create different 
wind and temperature regimes during various parts of the 
Gay. Just after sunrise, the heat received at the ground is 
transferred upwards in many ways. In a very thin layer above 
the ground, heat 1s transferred by molecular conduction 
along a very large temperature gradient. Above this, heat is 
transported almost entirely by forced convection. At even 
larger distances, a free-convection layer begins to form in 
which buoyancy production dominates. In order for the tem- 
perature of the air to rise, the heat flux (which is 
directed upwards so that 6@’w’ > 0) must decrease with 
height. According to Businger (1973), the heat flux at the 
Surface and the height at which the flux is zero generally 
increase until midday, ultimately resulting in a constant 
flux layer up to 100 m or so. Much later in the day, as the 
ground begins to cool through radiational losses, a downward 


heat flux exists and a stable layer forms close to the 
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ground. This represents the time of interest in this study. 
Near the ground, the flux Richardson number is usually small 
due to the large shear in the wind, and increases from the 
surface up to some maximum. In the upper part of the boun- 
dary layer, instability persists, resulting in a small or 
negative R;. The critical R, is reached first at the height 
of the maximum, and there a laminar layer forms, effectively 
preventing further exchange of heat or momentum across this 
layer. Temperatures near the ground decrease even faster as 
the downward heat flux is eliminated. Friction with the 
ground acts to slow down the air below the laminar layer and 
the associated reduction in shear causes R, to exceed its 
critical value, suppressing turbulence everywhere. Vertical 
transfer drops to near zero as a result. 

The time of interest in this study is after sunset, 
when the inversion is present throughout the valley and 
Slope winds are beginning to develop. For the reasons 
outlined previously, the vertical fluxes of heat and 
momentum are quite small during this period. As a result, 
the vertical transfer coefficients were set to zero in model 
integrations attempting to simulate realistic valley 
conditions. On occasion, the model is executed using neutral 
Stability conditions initially, and, in this case, the ver- 
tical transfer should be non-zero. However, definition of a 
vertical transfer coefficient which is a function of 
Stability was precluded in this study due to lack of time 


and lack of proper knowledge of its variation within small 
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valleys. This is left for a later study. 


2.3 The Assumption of Incompressibility 

The assumption of incompressibility of the flow is made 
(dp/dt = 0), resulting in simplification of the equation of 
CONELOULEY Faia. 22s) te AGCOLGINds tos Batcheloryig 967), ethe 
following three conditions are necesSary in order for air to 


behave as if incompressible: 


U2/c2 << 1 : C2. Gens 
n*L2/c2 << 1 ‘ G2is37 2') 
pgL/yp << 1 ; (Zeon) 
where c = speed of sound, n = meaSure of the dominant 


frequency of oscillations in the flow field, y = Cp/Cy- 
Spatial variations of V are characterized by some length 
scale L and variations in the magnitude of V with respect to 
position and time have the magnitude, U. 

The first equation, (2.3.1), states that the flow vel- 
ocity must be small compared to the speed of sound. For air 
at 15 °C and 1 atmosphere pressure, c = 340.6 ms~', and 
therefore (2.3.1) was satisfied for all flows considered in 
this study (where U is about 1m s~'). The second equation, 
(2.3% 2)omStatess thatathe frequency of any. oscillatory motion 
in the air must be small compared to the frequency of sound. 
The highest frequency resolvable in a finite-difference grid 
model is n = 1/2At and, in the cases considered, At was of 


the order of 1 second. Thus, (2.3.2) was always satisfied. 
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Equation (2.3.3) restricts vertical scales of motion of 
the air to much less than the *scale* hei ont! ,fHy of the*at= 
mosphere, which is the vertical distance over which pressure 
and density diminish by a factor of 1/e, and is equal to 
p/eg. For air at 0 °C, H = 8.0 km. Even when the temperature 
is not uniform, it is evident that this condition will be 
satisfied aoe Be motions occurring in layers of the atmos- 
phere not exceeding a few hundred meters in depth. From the 
application of Batchelor's criteria, it can be seen that the 
assumption of incompressibility is reasonable. Equation 
(2.2.3) can now be rewritten as 


VeV = 0 : (e354) 


2.4 The Vorticity Equation 

As 1S customary, the vector equation of motion is 
replaced by a scalar equation, known as the vorticity 
equation, by taking the curl of (2.2.1), or equivalently by 
taking appropriate Spatial derivatives of the two scalar 
component equations of motion. The latter derivation is 
presented here. The component equations of motion for v 
(cross-valley wind component in y-direction) and w (vertical 
wind component in z-direction) can be written as the 
following, after the previously mentioned assumptions have 


been» made® tow'( 232259 >% 
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ow ow Ow 1 9p 92w , o2w (2.4.2) 
AV + w s- + K, e+ of. 


where Ky = horizontal diffusion coefficient, and K, = verti- 
cal diffusion coefficient. Differentiating (2.4.2) with 


respect to y, and (2.4.1) with respect to z, and subtracting 


yields 
26 DE, 9E, - (av, aw} _1 (0 ap , a9 ap a2 cent ai ase 
yy By y My fy Aras es SE ee 
where 
ow av 
Ee ee Zot is 
Aare” ae ( 3) 


—& being the vorticity component along the down-valley axis 
(— < 0 for clockwise flow and vice versa). Using the 
continuity equation approximation for an incompressible 


Plusd, (2.3.4), one obtains 


CA tee pl sayy te Pad Sa Sh aged LS ae ae 
Yo aa Ce Zz 1, [fe 32 + 22 oy 1 Le cena dz (2.4.4) 


Therefore, in this model, the local change in vorticity, 
must be a balance between advective changes in vorticity, 
the so-called 'solenoid term' given by 
Wile OBAOV40D/AZmrn0PL0ShOR/ Ov peanGpGittusionoLavonticity. 
At this point, the solenoid term will be discussed in 
terms of how it is able to alter the vorticity values. Such 
a consideration may shed some light on the kinds of airflow 
or circulation expected to occur within the valley. In a 
barotropic atmosphere, p = p(p only), and there will be no 
contribution to (2.4.4) through this term. This can be seen 
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as given below 
C= fv dy + w dz 


where € = C per unit area. It then follows that 


dc dv dw 
dt dea ari cage ss 


If only the terms involving pressure gradients are retained 


in the expressions for dv/dt and dw/dt, one obtains 


dC i 
mm $= 


dt p {dy 


[BB ay - Bac] = - gd ay (2.4.5) 
With p = p(p), this becomes the closed line integral of an 
exact differential and as such is Zero. 

In general, density 1S not determined simply from the 
pressure (in which case the atmosphere is called baroclinic) 
and the isobaric surfaces (dp = 0) and the isosteric sur- 
faces (dp = 0) intersect each other. Consider the baroclinic 


fluid shown in Figure 2.1. Evaluating (2.4.5) around the 


path indicated results in 
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This results in growth in circulation in the negative 
(clockwise) sense about the curve, i.e. in such a direction 
as to let denser fluid sink while lighter fluid rises. 

The direction of the circulation which grows from a 
State of rest aS a nesuit of) bancclinicity turns the 
isosteres more nearly parallel to the isobars by moving the 
dense fluid towards high pressure and the lighter fluid 


towards low pressure. This process acts to convert potential 
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Figure 2.1 Relation between isobars and isosteres ina 
baroclinic fiuid 
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energy of the mass distribution of the fluid into kinetic 
energy of circulation, thereby minimizing the potential 
energy. 

In this model of the boundary layer, as in many others, 
e.g. Malkus and Witt (1959), Thyer (1966), and Orville 
(1964), the solenoid term was approximated to be (g/@ 06/dy) 
The equation of state, (2.2.4), by taking logarithms and 


differentiating, can be rewritten as 


L opeenl epee 1 ot 
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These are used to rewrite the solenoid term as 
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Adding and subtracting g/@ 06/day with 


from the Poisson equation, (2.2.5), and simplifying, results 


in 
slept le epee lf er | tuel Re lop: 
pT az dy ' pT ay az 6 dy & {T ay cP by 
; aT 
Letting Y = - 55 = ambient lapse rate C26) 
and Yq See dry adiabatic lapse rate (2747) 


P 
one obtains, after combining terms, that the solenoid term 


can be expressed as 
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(A) (8) (C) 


(2.476) 


Rough estimates of the magnitudes of the terms in 
(2.4.8) can be made by considering the horizontal and verti- 
cal pressure gradients which exist within a river valley 
such as the one being considered in this study. Up to this 
point, no mention has been made of the hydrostatic 


assumption, 


8Ph P 
Pe ee ee a4. 
dz RT ! = 


obtained from the third equation of motion, (2.4.2), when 
vertical accelerations and diffusion are neglected. The 
question of validity of this approximation can be found in 
comparing the magnitudes of typical hydrostatic and 
non-hydrostatic pressures. 

First, hydrostatic pressures are computed at two points 
within the valley where differences are expected to be 
largest. The hydrostatic equation above, (2.4.9), can be 
used to integrate downwards to estimate the pressures 
resulting from various temperature regimes in different 


parts of the valley using 
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The assumption is made that at z = z’, the temperature, T 


’ 


and pressure, p , are constant. 
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Based on observations within the North Saskatchewan 
River valley on 21-22 October, 1977, horizontal temperature 
differences between the slope and the free air of close to 
1 C° were recorded over a distance of 185 m (Paterson 
(1978)). These values were obtained from a thermograph 
located on the slope 26 m above the valley floor (station 
#2) and from averages of thermocouple readings obtained in 
the free air at heights of about 30 m (station #4b) and 22 m 
(station #4c) above ground. The latter were suspended from a 
bridge which spans the valley. The locations of the 
instruments are shown in Figure 2.2. Some assumed tempera- 
tures and corresponding lapse rates (which were constant 
with height) were used to simulate the above observations. 
The pressure and temperature at the top of the model 
(z° = 108 m) were assumed to be 917.8 mb and 278.632 °K, 
respectively, as assumed in the standard run - see Chapter 
4. The variation in temperature with height is given by T(z) 
= T° + 7(108 - z(m)). Figure 2.3 illustrates these 
assumptions. At location A, (2.4.10) becomes 


BE pk <i ga(to digg Beles deh 
P R R T’ + y (108-2) 


Pa 26 26 
or 
Tee ln T Ie 
= ra SEEM Ao IST 
Pa T’ + y (82) 


where the integration on the right hand side has been 
performed using a change of variable technique. Evaluating 
the above, using a horizontal temperature gradient of 1 C° 


per 200 m, one obtains xk 927.097 mb and Saas 927 ..079 «mb: 
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Figure 2.2 Location of instruments used to estimate horizon- 


tal temperature differences within the North Saskatchewan 
River valley 
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Figure 2.3 Assumed temperatures and lapse rates used 
estimate hydrostatic pressures at A and B 
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This results in a horizontal hydrostatic pressure gradient 
Ofe79) x Os °¢ mb mii'r OTe Sark 10g *) Pa- ms)... Repeating this 
calculation for a horizontal temperature difference of 
10 C°, where Tac 267.992 ‘Kegand Ty = 277.972 PKs Yn 
=a 0 13510; Gre Nea ena ae =O; 01084 Cro Meng HESUd tS; Wlinsae Nonizon= 
tal hydrostatic pressure gradient of 6 x 10°? Pa m-'. This 
lis expected to be an overestimate of the real gradients 
involved. 

Secondly, the non-hydrostatic pressures at these same 
two points are considered. Equation (2.4.2) can be 


approximated as 


neglecting diffusion. By consideration of the typical 
accelerations involved at the positions under consideration, 
A and B in Figure 2.3, estimates of Py, can be made. Close to 
the slope, downslope winds up to 0.5 ms~' are typical in 
the North Saskatchewan River valley. Assuming the angle of 
the sloping ground to be 0.1 radians, a vertical velocity 
component of about 0.05 m s~' results. Thus, w undergoes 
some deceleration from 0.05 ms~' at mid-slope (say, z = 25 
m) to 0.0 at the horizontal valley floor. This occurs over a 
time interval equal to about 500 seconds (the downslope 
wind, having speeds of 0.5 ms~', covers 250 m - mid-slope 
to-syanl l exact Looteury ein: abouts 500..seconds.)). Substututing this 
estimate of dw/dt, and integrating dw/dt = -1/p OD ,/ Oz from 


Zea. tO? Zs 25 m, one-obtains 
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Pohl, = = dw/dt |, Az =i+ge xe 10 mb. Simrlar casculations 
within the free air, assuming w accelerates from 0.0 at the 
Vauley floor to ebout 0.005 m- si at’ z ="25 m in5000 sec-— 
onds time, result in Pnhla= +3 x 10°’ mb. These rough 
calculations reveal that the non-hydrostatic pressures are 
much smaller than the hydrostatic ones, and, as a result, 
the pressures involved in the study are clearly principally 
hydrostatic in nature. The horizontal gradient in non-hydro- 
Static pressure, which may contribute to horizontal 
accelerations through equation (2.4.1), are found to be very 
small as well. 

Returning now to the evaluation of the terms comprising 
the solenoid term, given by (2.4.8), one finds that term B 
is identically zero if the pressure is assumed hydrostatic. 
Terms A and C can be estimated based on the rough 
calculations given above. Term C depends on the horizontal 
pressure gradient as well as the stability of the air. Table 
2.1 indicates the relative magnitudes of these terms for the 
two cases considered above, assuming p = P,? mano rizon= 
tadm temperatures differencesofi 1 Co overs'200. m;, 2;)) axhord- 
zontal temperature difference of 10 C° over 200 m. Potential 
temperatures required for term A are computed from the 
assumed temperatures and hydrostatically derived pressures. 
From the Table, it is clear that term C is much smaller than 
term A, and therefore is neglected in the model. 

Pivnialily), theevorticity equation, (2.74.4), can be 


expressed as 
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Table 2.1 Estimates of magnitudes of term A (g/@ 06/dy) and 
term C (ap/dy R/p g/C, (1-7/y,4)) comprising the solenoid 
term 


aT/ay = 1C°/200m 
neutral (y=+.01 C° /m) 
stable (y=-.004 c°/m) 
very stable (y=-.04 c°/m) 


3T/dy = 10C°/200m 


neutral (y=+.01 c°/m) 


stable (y=-.004 c°/m) 


very stable (y=-.04 c°/m) 
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2.5 Use of Streamfunctions and Related Boundary Conditions 
For two-dimensional incompressible flow fields, the 

non-divergent velocity field, Vey : + Ww k, where 5 and k 

are unit vectors in the y and z directions, respectively, 


can be expressed using a streamfunction, ¥. In this case, 


the velocity, streamfunction and vorticity fields are 


related by 
v= ov (255.1) 
oz 
Vio! (25.2) 
dy 
and 


po - vy = -[F¥ + Hy) (2.5.3) 
These three diagnostic equations were used to obtain the 
velocity components, v and w, from the vorticity values 
predicted by (2.4.11). 

In order to satisfy certain velocity boundary condi- 
tions, corresponding boundary conditions must be imposed on 
the streamfunction, yw. The relationship between velocities 
and streamfunction may be better understood by writing V as 
-ix Tue where i ifatlunits vectorPanetherxeadirection=. it Vw 
happens to be in the y-direction, then V is in the 
z-direction, and vice versa. In other words, Vis parallel 
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Since axial symmetry was assumed (about the vertical 
lines bordering the valley cross-section briefly described 
previously), the horizontal velocity component was required 
to be zero. Equivalently, this disallowed outflow or inflow 
along these lateral boundaries. Therefore, dW/daz = 0 or Ww is 
constant with height there. 

Using similar arguments for the upper and lower 
boundaries, W was required to be a constant everywhere along 
the cross-section boundary. Since W occurred only in the 
form of derivatives in the equations, this constant was 
arbitrary and set to zero for simplicity. Flow everywhere 
parallel to the boundary ('free-slip') was guaranteed 


through these boundary conditions on the streamfunction. 


2.6 Potential Temperature Boundary Conditions at the Ground 
The diurnal temperature change at the surface is the 
driving force for the valley wind development. Cooling along 
the valley slopes creates a non-zero horizontal potential 
temperature gradient, and this in turn is responsible for 
the creation of vorticity (through the solenoid term). 
Following the approach of Thyer (1966), the surface 
temperature boundary condition was initially determined 
through a requirement that the temperature gradient at the 
ground surface be consistent with the rate of heat flow and 
the eddy conductivity. Assuming some rate of heat loss per 
unit area, Qo, the outgoing energy from the sloping surface 


is Q, cos (a), where a is the angle the sloping ground makes 
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with the horizontal. In order for this heat loss to be 
completely balanced by conduction from the air with some 


eddy conductivity, K a certain vertical temperature gradi- 


Te 
ent must exist which satisfies 


08 
eae Me CTIO page (2.6.1) 


based on the heat conduction equation. 

Such constraints, however, provided inappropriate cool- 
ing rates at the surface. In reality, surface temperatures 
are altered through radiation effects, advection, and diffu- 
Sion processes, and this early approach was soon abandoned 
in favour of simply prescribing a rate of surface tempera- 
ture change. 

Some observations of surface temperature changes within 
valleys are available. Orville (1964) chose to model the 
diurnal trend in temperature in mountainous terrain as 
reported in Geiger (1957). In his model, a Sinusoidal pot- 
ential temperature change at the surface was used, having a 
period of 24 hours and an amplitude of 7 C° at the valley 
bottom and decreasing linearly with height to 3 C° at the 
top (1 km above). 

Rao and Snodgrass (1981) attempted to determine a 
steady-state drainage flow, and they used a constant cooling 
rate of 2 C° per hour for the first hour. Subsequently, the 
surface temperature was unchanging. Steady drainage flow was 
attained "in a few hours". The specified cooling rate was 


rather arbitrary, and the time scale of the evolution of the 
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flow was interpreted simply as an indicator of the 
corresponding surface temperature deficit. 

In most valley model integrations in this study, which 
lasted well under an hour, a constant Surface cooling rate 
of 2 C° per hour was used. All the sensitivity analyses were 
performed under these conditions. A radiative approach to 
the problem of heat transfer in the lower layers of the at- 
mosphere was used by Brunt (1934). He showed that radiative 
transfer of heat could be approximated as a diffusion 
process, with a constant exchange coefficient, Kp, called 


the radiative diffusivity, i.e. 


a 9°T 
ne SR ae C261, 22) 


This is based on the absorption spectrum of water vapour ; 
long-wave radiation emitted by the surface of the earth is 
partly absorbed by the water vapour in the atmosphere and 
reradiated. The numerical value of K, was a function of tem- 
perature, the distribution of vapour pressure with height, 
etc. but by using typical values, Brunt suggested that 

Kp = 0.13 m* s~', Anfossi et al. (1976) experimentally 


detrerminecsa:valuesof 0.3:me si” 


under conditions favourable 
for radiative inversions. The use of a diffusion process to 
model the radiative heat transfer near the ground was 
attractive and easily implemented into the model under 
development simply by modifying (K,), in (2.2.2). 


Ditferentiating (2.6.2) once witherespect. to 2) .and 


replacing oT/dz by the thermal flux, f, one obtains 
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afl, tt (Ges) 
The solution to this differential equation was obtained by 
Brunt in order to predict nighttime minimum temperatures, 
resulting from radiative cooling at the ground alone with a 
eonstantaneat flUx;eio: <8 01/0z),—5. In other words, it was 
assumed that the net loss of heat by radiation from the 


ground to the atmosphere, R,, was constant. From (2.6.2), R 


N’ N 
can be expressed as Kp p Cy, 0T/0z|z=20- By integrating the 
SSLUMI ON=toB C256 SS) FSECZ ete eMthe original “variable; *T( 27) , 
was determined to be 

= Un 1s Zee 
a2.) el 0. 0) = — (K,t) “exp (-27/4K, t) - | exp (-u2) du 
noC K. Zz : 
P hy CZ obi 4} 
2(K,t) 
Therefore, the surface temperature, T(0,t), is expressed 
Simply as 
5 
2 Ry t 
TRO) as TO rararelan (2.6.5) 


25 
C 
(1K, ) “p 5 
In Chapter 4, the results of an integration using (2.6.3) as 


the surface temperature boundary condition are given. 


257| Initial Conditions 

For the initial conditions of a state of complete rest, 
the vorticity, velocity components, and streamfunction were 
all set to zero. The only other quantity which required ini- 
tialization was the potential temperature. Given an initial 


lapse rate and the temperature and pressure at the lowest 
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level, the potential temperature field everywhere in the 
‘valley cross-section was computed using Poisson's equation, 
(2.2.5), the hydrostatic equation, (2.4.9), and the equation 
of state, (2.2.4). Poisson's equation can be written in 
logarithmic form as 


d 1n6é = d InT -2 d lnp 
P 


Substituting (2.2.4) and simplifying results in 


dz }dT . B 
d iné 7 E + i] 


Using the lapse rates defined by (2.4.6) and (2.4.7), one 


obtains 


dz 


d ind = > fy, - y} (22 Teas) 


The potential temperature at the valley trough can be 
computed using (2.2.5) given the trough temperature and 
pressure. A surface pressure of 930 mb was used for 
Edmonton. The potential temperature at any height z above 
the valley floor is computed using (2.7.1) and by 
integrating upwards from the ground surface. The actual tem- 
perature, T, is required for this computation but it can be 
found from the known constant lapse rate. The isentropic 
(constant potential temperature) surfaces are assumed 


constant with height initially everywhere within the valley. 
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2.8 Summary - Prognostic and Diagnostic Equations 


The simplified prediction equations which describe the 


thermodynamics and circulation are 


96 00 96 9269 926 
ae Vay “az * By yz + By gaz ee 
dt dy 2a Ody: y oy* z oz (2.4.11) 


Given initial and boundary conditions, these equations were 
solved numerically. From the & field so derived, the 
following diagnostic equations were used to solve for the 


corresponding velocity components : 


ay 
“its 627553 ) 
ORCA 
w = dy (245°.1) 


E= - Vy (25522) 
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Chapter 3 


NUMERICAL ASPECTS OF THE COMPUTER MODEL 


3.1 Introduction 

As previously indicated, the non-linear character of 
the system of equations to be solved necessitates the use of 
numerical methods, since general analytical solutions are 
not available. To this end, the relevant equations were set 
up in finite-difference representations in time and space. 
The numerical method described is, in principle, applicable 
to three space dimensions as well as to two. However, the 
actual computations are two-dimensional in order to make 
computer storage and run-time costs feasible. According to 
the Eulerian system of equations, time (t) and two space 
co-ordinates (y and z) were chosen as independent variables. 
Boundary conditions appropriate to the chosen finite-diff- 
erence formulation of the equations are required. 
Approximate finite-difference formulae to the equations 
under consideration may be obtained by various methods. 
However, each formulation possesses certain levels of 
accuracy, consistency, stability, convergence and 
efficiency. 

The accuracy of a given finite-difference formulation 
is determined by estimating the truncation error involved in 
estimating a aaa ae Say, by some grid point 
approximation. For such an approximation to be considered 


consistent, it should approach the derivative as the grid 
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interval approaches zero. Given that the true solution is 
bounded, a finite-difference formulation is termed stable if 
the difference between the numerical solution and the true 
solution remains bounded with time. Typically, some 
constraint on the time step, At, 1S required, in order for 
the method to remain stable. If the error approaches zero as 
both the grid size and time step are decreased, then the 
solution is called convergent. For properly posed initial 
boundary value problems described by partial differential 
equations with consistent finite-difference formulations, 
Stability is the necessary and sufficient condition for 
convergence. Finally, computer storage and costs must be 
considered in the selection of a numerical scheme. For a 
Simple problem with a small number of variables, little 
computer time is required and a Sophisticated scheme of 
high-order accuracy may be used, whereas for a system of 
many variables, accuracy must be compromised in the 
interests of efficiency. Whether a given scheme is explicit 
Or implicit comes into consideration in this regard. For an 
excellent discussion of numerical methods for atmospheric 
studies, see Mesinger and Arakawa (1976) or for a more 
mathematical discussion see Potter (1973). 

In order to compare the numerical solutions obtained 
from one ate finite-difference formulations for diffusion 
type terms, a simple equation for which the analytical 
solution was available, was used. This was the one-dimen- 


Sional heat conduction equation, subject to certain relevant 


a 


™m 


eit noituloe sutd om 


ai nots 


*p -— = 


=a 


— 
_ y 


noisulos teotysqu 
vilsstqy? ake dtiw & 


“sg 2a1re af9.41 sidsta the 


S16 get2 emis Bane 


L 
| saidese Semzes 
’ 
sut? of3 Bris 
STO 
sriupes er 
i a —~/; 
‘ é \; few 
ad 592437 278R 
a“ + CT ¢ I 
| “te 
Pp : 
} 
: 
| 7 I { 
:¢ 
, 
‘ - 
rs 
j T= 
rf ( 
i% f 


54 t2etdve \nebeeupe pee ape: £4) 


2 


fr 


th .qese ae 


t. 
— 
as 


“a 


eb ametdone. adil 


es 


Sasven ot 


estar ‘sneteras 


cos 


~~ e _ 


Ctttd <.8He 6 skeen 


= 
5 BAS Soi tupes ai og 
_ 


Se ai yer YoOsuvsoe: 


wo 


ntyvoDe oa 


= 


arehiena> oank conta sal 
soivemun 10 notesvosshu 
SABA OAs 336 

trod sea noizevoet 

nin sdf e2agmoo O47 1a670 4 

: sonesetiS-setiati eugssa 


cigmia s.,2 


Pics 


ldsiteve. egy 


ri > 


_ 


is 


—~ 


38 


initial and boundary conditions. In this chapter, various 
numerical solutions to this equation are presented along 
with the analytical ones derived. In the early stages of the 
valley model integration, vertical diffusion processes 
Gominate, and, as a result, the solutions obtained from the 
one-dimensional diffusion equation with the same initial and 
boundary conditions are very similar to those of the valley 
model equations. This provides a check on the early stages 
of the evolution of the thermodynamic fields. Finally, in 
this chapter, a technique for solving an elliptic Poisson 
equation is outlined. It was used to obtain the streamfunc- 
tion field, given vorticity values, through equation (2.5.3) 


and is termed Liebmann sequential relaxation. 


3.2 The Grid 

In order to study the circulation within a small urban 
river valley, it is necessary to define a physical boundary 
which is realistic but also fairly simple. A nearly 
'V-shaped' valley was used with physical dimensions 
resembling those of the North Saskatchewan River valley. Due 
to the meandering nature of the river, asymmetries in the 
topography of the valley are evident. Figure 3.1 indicates a 
typical vertical cross-section through the valley, close to 
the city center, at a point where the river is aligned 
east-west. The north-facing, or south, slope was modelled 
in this study, using a much simplified lower boundary 


indicated by the dashed line. 
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Most simulations, therefore, occupied a two-dimensional 
region bounded by a lower surface consisting of a ridge 52 m 
high, and a horizontal valley floor 80 m wide. The slope of 
the ground surface was approximated to be 0.1, forming an 
angle, a, of about 5.7° with the horizontal valley floor. 
Some results are also presented in Chapter 4 for a slope of 
0.2. In order to associate grid points at the boundary with 
the surface of the earth, constraints are placed upon the 
horizontal and vertical grid spacings permissible in that 


region, Ay and Az, respectively such that 


oan @) = = (Baza) 


Close to the ground surface, gradients in most 
variables are large and higher resolution in terms of the 
grid spacing is needed. A logarithmic grid was felt 
difficult to implement due to the geometry of the area being 
considered. A finite element approach, which can incorporate 
variable grid spacings and is generally well suited for 
complex terrain, is now becoming useful in meteorological 
applications. However, it was not considered at the time 
when model development was begun. Use of terrain-following 
co-ordinates has much to offer in terms of reducing complex 
geometry to a simple rectangular region but introduces many 
other complications once the equations have been 
transformed. AS indicated earlier, a finite-difference 
technique was selected to integrate the equations numeri- 


cally. The grid length finally chosen was constant in the 
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lower part of the valley cross-section (Az,) and larger, but 
still constant in the upper region (Az,) where the variables 
are expected to change more smoothly. A constant grid length 
(Ay) was used in the horizontal, or north-south, direction. 

According to Thyer (1966), the height of the horizontal 
upper boundary should be at least twice the ridge height, if 
Dees Ato include the entire region which such a circulation 
is likely to occupy. In most runs, the height of the upper 
boundary was chosen as 108 m but sensitivity to this value 
is considered in Chapter 4. A schematic diagram Of the grid 
DoOimmtlattice 1S Giveniwin Figures 3.2. 

At any time t, any point in the region can be referred 
to by a set of indices (J,K,n), the points being numbered 
horizontally by J from the valley axis (J = 1) to the ridge 
line (J = JMAX), vertically by K from the valley trough 
(K = 1) upward to K = KMAX, and a time index n where 
t =n At. This notation will be used in subsequent sections 
which deal with the finite-difference formulations and boun- 


dary conditions. 


3.3 Initial and Boundary Conditions in Terms of Finite 


Differences 


3.3.1) Initial Conditions 
Since the atmosphere is assumed to be at rest ini- 
tiatly, allyorrcepoine values ob vorticity, §. streamfunc- 


tion, Ww, and velocity components, v and w, are set to zero 
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Figure 3.2 Grid configuration for valley model 
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E(J,K,0) = 0 


¥(J,K,0) = 0 l1<J< JMAX, 


(Sarserih) 
v(J,K,0) 


| 
(o) 
~~ 
‘N 
w 
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w(J,K,0) 


a 
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The potential temperature at the valley trough, 6(1,1,0), is 
easily obtained from the initial temperature, T(1,1,0), and 
pressure, Dil, l,0); using (2.2.5), tor be 


R/C 


6(1,1,0) = | P 7(1,1,0) f eee 


by 
caictaeh 200 
PC ii3 0) 


Remember T(1,K,0) is specified by the valley trough tempera- 
ture, T(1,1,0), and the known lapse rate, y, through 
(2.4.6). Using (2.7.1) to obtain the potential temperature 
field directly above the valley trough by integrating 
upwards uSing centered differences, one obtains 


6(1,K,0) = exp{ln 6(1,K-2,0) + Cy, 7 Y) 2 42/T(1,K-1,0)} ES ASHE ASS) 


For the first interval, a one-sided difference formulation 
must be used, using the mean temperature for that layer, 
CT (1 we 0.). tecT fal. 2. Ow) AZ, nas, below, 
6(1,2,0) = exp{ln 6(1,1,0) + (Yq - y) Az/ 
one eee) 
(Tl, 1,0) et VY 5250) 2) 
The isentropicesurfacesyare all horizontal initially, and, 
therefore, 


6(J,K,0) = 6(J-1,K,0), l1<J<«< JMAX : (530098055) 


and the entire potential temperature field is now specified. 
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3.3.2 Lateral Boundary Conditions 
According to arguments outlined in Section 5 of Chapter 


z, the flow must be purely vertical at the lateral 


boundaries, J = 1 and J = JMAX. Therefore, 


v(1,K,n) = 0 1 < K < KMAX, 


(3434.6) 
v(JMAX,K,n) = 0 n 20. 


As well, these boundaries form lines of symmetry, and the 
circulation on either side is assumed to be a mirror image 
of that computed within the valley cross-section being 
considered. Thus, the streamfunction values on one side of a 
lateral boundary are equal in magnitude but opposite in sign 
to those on the other. See Figure 3.3. The vertical velocity 
along these boundaries, e.g. at the grid location denoted by 
BeChmEGUrC 2 wey Canebe Computed irom (205° 2) assuming 
Wis) p= = yt); to be 

w(A) = - (CC) - p(B))/2dy = - p(C)/dy 


In general, then, 
w(1,K,n) = - y(2,K,n)/Ay 1 < K < KMAX, 


w(JMAX,K,n) = p(JMAX-1,K,n)/Ay ) n 2 0. paras 


Since dw/dy = 0 through symmetry and v = 0 as before, 
the vorticity along the lateral boundaries must also be 
Zero sthrough. (25493) 2) Thus; 

E(1,K,n) = 0 1 < K < KMAX, 
| (3.3.8) 
E(JMAX,K,n) = 0 ) n 20. 


The streamfunction is constant and arbitrarily set to zero. 


y(1,K,n) = 0 1 < Ks KMAX, (3.3.9) 


w(JMAX,K,n) = 0 ) n 20. 
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Figure 3.3 Schematic illustration of axial symmetry and 
boundary conditions 
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Horizontal temperature gradients across the lateral 


boundaries vanish through symmetry. 


3.3.3 Lower Boundary Conditions 

As free-slip boundary conditions are applied at the 
ground, velocities must be parallel to the lower surface. 
Along the horizontal part, therefore, w(J,K,n) = 0, and v is 
determined from the vertical streamfunction gradient using 
one-sided differences, since streamfunctions below the 
ground are not available. A one-sided formula, using two 
points above the surface, is given by 
v(J,K,n) =o 


eZ (5,50) (Beas 10) 


= {-3/2y (J,K,n)+2y (J, KF1,n)-1/2¢(J,K+2,n) }/Az 


This approximation is also used to obtain the vorticity, 
Since & = - dv/dz there, or 
E(J,K,n) = -{-3/2v(J,K,n)+2v(J,K+1,n)-1/2v(J,K+2,n) }/Az ee eoe: Hee) 
Along the sloping ground surface, the magnitude of the 
velocity is computed from the streamfunction gradients above 


the slope (again, in a one-sided fashion) since from (2.5.4) 
ss ar 1 
[V| = [Vo] = {@ep/ay)? + (ay/az)2}* (oear tcp 


The individual velocity components are then determined using 


Vevjt+twk=-(|¥| cos(a) 3 + |V| sin(a) k), where 


R 
i} 


angle of the sloping terrain. The minus sign appears in 
order to ensure downslope flow. In finite-difference form, 


then, 
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Al = {-3/2y(J,K,n)+2p(J-1,K,n)-1/2yp(J-2,K,n) }/Ay : 
visa) 
oe | = {-3/2p(J,K,n)+2y(J,K+1,n)-1/2y(J,K+2,n) }/Az , 
Zz 

(34K,0) 


v(J,K,n) =-{ (80/291 05 x ny)? + (BV/8z| (5 ny )?}? cos (a) é 
CBsSenlise) 


w(J,K,n) =-{ (av/ay| (, K Lr + (3y/82| (; ra aor sin (a) 


See Figure 3.4. The velocity components, at grid points 


which intersect the sloping surface, must satisfy 


wee = tan (a) = slope of ground (E3202) 


and flow parallel to the material surface is thereby 
Quaranteed. Having obtained the surface velocity components, 
the surface vorticity is evaluated using (2.4.3), from 
one-sided differences of the velocity components, as usual 
E(J,K,n) = {-3/2w(J,K,n)+2w(J-1,K,n)—-1/2w(J-2,K,n) }/ay 
(Sig She Sy, 
- {-3/2v(J,K,n)+2v(J,Kt+1,n)-1/2v(J,K+2,n) }/Az 
The streamfunction is set to zero on grid points 
intersecting the lower boundary. The potential temperature 
is Simply prescribed as some function of time, as alluded to 


earlier. 


3.3.4 Upper Boundary Conditions 
Consistent with the upper boundary being located above 


tne valley circulation, both velocity components vanish 
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Figure 3.4 Schematic illustration of lower boundary 
condition 
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along K = KMAX. In addition, the vorticity and streamfunc- 
tion values are set to zero. The boundary conditions on the 
flow are Summarized in Figure 3.5. The potential temperature 
at the top of the region is assumed constant, being far 


enough away from the cooling surface. 


3.4 Finite Difference Formulation of Prediction Equations 


3.4.1 General 


The simplified prediction equations to be formulated 


using finite-differences are 


90 = 38 36 926 370 
ate 4Y dy he ap i pees ree $ (2.2.2) 
Seay ewe eK are eK "5 4 2 38 (2.4.11) 
at ay moze) yay zdz- 6 dy 2.4.11 
advection diffusion other 
(a) (b) (c) 


Typically, it is possible to examine the stability 
properties of a given numerical scheme only when it is 
applied to simple types of differential equations. This is 
because it is only with simple equations that the 
mathematics involved become tractable. When an equation 
involves more than one type of term (e.g. advection and dif- 
fusion, as in (2.2.2) and (2.4.11)) Mesinger and. Arakawa 
(1976) recommended the use of "different schemes for the 
different types of terms". For this reason, terms (a), (b), 


and (c) were considered separately and the most appropriate 


oS ae 
“sauimesti2 bas yti ; 
sqi no enoizibnes ye 


fettasteq oa? . 8 omy 


4 
] 
| 
i 
- 


. 4 
= 
: 
a 
S 


- - 3 . 7 - a % 
noisazibertd %% nolisiumtetl soaeTte7zeze 
_ 


_ 


‘sups nersor pas bs its igaks 
516 enscigaeti rae 
es ; . 


“+ 
i+ oS WV Se ee 
° 


2 


ame apivemun nevip 210 

feisas Tei tis Be zaqy? siquia 03 
itewps eigmie dsiw yine er 

,22 amoned beviognt- 
y¢ 20 eagys sao nals Sree 

e > 
12f290M (érp .By8) Hite (¢.$.¢) ak #6 
4 
bib” Io. sau end arereen 7 


sm1es ,fo8er :iay 40%. "smi89 ic sole fete 


1% Jzom oad one iesssegee bereb ia: 


a ae 


vew= Eun yp = 0 


Elen ima ims 


w = p(JMAX-1,K)/Ay 


Vee cen yeai0 ; 


w= -(2,K)/Ay 


< 
n 


v= ay/dz, * (v| cos(a), * 


ka sin(a), * 


fe 


w= y = 0, 


aw/dy - dv/oz, * 


wr 
L] 


—E = -dv/dz *, 


| 


y= 0 


* one-sided differences 


Figure 3.5 Summary of boundary conditions on velocity 
components, v and w, vorticity, &, and streamfunction, yw 
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numerical scheme chosen for each. As it is difficult to 
determine the stability criteria for non-linear equations, 
it was assumed that linear stability criteria would provide 
a close approximation when applied to the prediction 


equations. 


3.4.2 Advection Terms (a) 
A comparison of numerical schemes available for solving 


a one-dimensional advection equation, 


——+ Co-— = OF c 230 $ (341) 


is given in Table 3.1. After examining the advantages and 
disadvantages of each, the forward-time, upstream-space 
formulation was chosen to determine the changes in vorticity 
and potential temperature through advection. Aside from its 
explicit two-level nature, this scheme was recognized to be 
reasonable in a physical sense, at least intuitively. 
Gradients in S upstream of the grid point under 
consideration are advected towards that grid point whereas 
those downstream do not enter into the computation. This 
method has been recommended by many authors such as Orville 
(1964, 1965, 1968), Molenkamp (1968), Pielke (1974), and 
Mesinger and Arakawa (1976). According to stability 


considerations, the solution will remain stable if 


A 
JEL Ee (3.4.2) 


< 


Ay 


However, damping proportional to wavelength can occur and is 


,enotteyps iussnii-aem 


shivorqg bluow si 29dPse 


os sluni2ikS af 


dan pes 


ot 


no f3i6 oq odd oo Baiigge © 


(s) ome 
ss 
>2 1072 sldalievse Saaenige fentaemua Fo As 


no bteups aoisoaves 
oso napa eee 


scatsevbs ad? gooinimexe wsttA <tek pid: 

ej sge~tasizaqu ,smie-biswi0d etd .aaee is. 
is10v ni asenreio oft snimesse6 oF AseoeD 
noit shied .sof?oevhs dovalld stuseTegees 
5+ besinpess: eew amedce aide ,ateten Deve 
vievisiueni seepl de .senee Deneayag & 7 
steoc btw ef7 to teeToagY : of ¥ 


dail 


he 


isdw jsniog Gss6 paces spuares S67 >evbs we not 
( 8] 
sidT .noisezuames ef? elant esis 3Cn ob MSO: 
re 


c 


rv10 es dove e1edtus ynem yd Betrenmtons7 
sng .(aCO)) sdfete , (286?) qwetactom , (6aRF 
vailtidese ot patbyenséA .(372!) svesers bas % 


sidate aiswe: [Liv noisuioe edd ,@nor: 


~~ fit! a 


vs eet 
-] bos w>:e nao doprelevew oF Eenodsicqemg gniguad \7 
} 


ey wt” : 


52 


uotzenba uotzdaape ue of patt{tdde SeweyoS Teostisunu swos jo satiiadoig L°€ aTqey 


(8961) dweyuaToW st yoTuM y ydeoxe *(QO861T) SWETTTIM pue AsUuTITeH : soUsIAazoy 


(squyod pyid Auew se 


SowT} 7) ot} Aaynduos pasearouy 


(stia} 


d10W) swt} Jayndwod pasearduUT 
*suba *3[NUTS 0} UOT ANTOS 

XV Z 3e Burdwep vwa1ijxe 
(a8ie[—°dwod) uoft4njtos uf sapom 7 


aut} Buy3ndwod peseeidutT 


‘yq3uaTeAem uo spuadap 3utduep 


AVETFqeqs ou - pastapeun | 


yasuaToARemM 


uo spusdsap 3utduep 


*dxa ‘ TaAaT-z 


> dxa  TaAeT-¢ 
"dur ‘ToAaT-Z 
*dxa ‘ [aAaT-Z 


°dxa ‘ TaAaT-¢ 


°dxo * [aAaT-z 
*dxa * TaAoT-z 


aFOFTTdxe 


‘[aAeT-Z 


| AVFTFQEAS | Anuayoyssa 


(,,4V+73V7)0 


(,,4V+z3V)0 
(,Av+23V)0 
(,Av+237)0 


(,Av+z3V)0 


(,4V+37)0 


(,Av+3y)0 


(Ay+37)0 


x SSTaM-Sj1eqoy 

aoeds 

lapio ya? SouwT}-— Ud 
Teptrozedelzy 
JJOIpUay—-xeyT 


aoeds—juo ‘ewtq—- uo 


(ouns je) 


piemyoeq—iatny 
aoeds—quo ‘owyt}-pay 


a3oeds 


| -weaizsdn Sout 3—pmy 


= 


PES Ss » 5 7 ~ F q22 3 


y — —_ ~~) 


- a. ~~ 14°. eo : 


- aalsl | "act 


7. z stables 


sidaseaxv | .gxo,Isyel-S 


= 
wy 


if gasioveaw ac shweqe’ an I , 2ats) ‘aie Sj (*obeaS96 | brewload-t4sfek \ 
: ve ; 
“sata gutteqne bonastont . . | is) 
i |: sh | ) 
> s: tasr6f-. quo) aoisufon at esi ; Mio! }. Gas is & {4° TA Jhj)0 pszege-2m , oa -~ $n 
ur ; i 
@ ; : ; 
; 53 S Ja a2atenns tia | 72094 msl f * 9} . tot -%1 , 


-enpe .2funla gnlaulec | arrel tye" 1490 3 Labiaasqet) 


evel eat? +9) i> ek | : ) | at-= 4 ‘oo byebrc | a 


Si 


most pronounced for wavelengths of 2 Ay when |c| At/Ay 
= 1/2. Further details can be found in Appendix A. The 
forward-time, upstream-space formulation, applied to a 


non-linear advection equation, such as 


of ae ot (37463) 


can be expressed as 


E(J,K,n+1) = E£(J,K,n) + At {A(J,K,n) + B(J,K,n)} 


where 

-v(J,K,n) EG. Ren) — GCL Rn) v(J,K,n) > 0 
A(J,K,n) = 

-v(J,K,n) (GH Kn) = $0, B0)) v(d.Ken) <10 

~w(3,k,n) ‘SCOK.n) — SOK) as,xsn) > 0 
B(J,K,n) = 

-w(J,K,n) {S(J,Kt1,n) - S(J,K,n)} w(J,K,n) < 0 


Az 


3.4.3 Diffusion Terms 
A comparison of numerical schemes available for solving 


a one-dimensional parabolic diffusion equation such as 
oS ix geome ny (34.4) 
dy ; 


is given in Table 3.2. Initially, the Dufort-Frankel scheme 
was chosen to evaluate the changes in 6 and £ due to hori- 
zontal and vertical diffusion. Even though this scheme 
requires computer storage of the vorticity and potential 


temperature matrices at three time levels, the extremely 
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Stable nature of the solution warranted its use. Further 
details can be found in Appendix B. The Dufort-Frankel 
scheme has been widely used in numerical modelling of the 
boundary layer (e.g. Estoque (1963) and most recently Kozo 
(1982)), and is favourably promoted in mathematical 
discussions of numerical techniques for solving parabolic 
partial differential equations (e.g. Richtmeyer and Morton 
(1957), Ames (1969), and Mitchell and Griffiths (i1980)). 

If a centered-time, centered-space formulation was 
applied to a one-dimensional diffusion equation, such as 
(3.4.4), it would be expressed as 

S(Jynr2)) ="S(J,n) ae S(irl ont lje= 2S (J.nti)) + S(J-4 nae) (22 Ash) 

Zt Ay ‘5 

where S(J,n) denotes the value of S at the Jth grid point at 
time t = nAt. This scheme would lead to numerically unstable 
solutions, as noted in Table 3.2. However, replacing 
S(J,n+1) on the right hand side of (3.4.5) by 
Usnue nat) SiO, n+2))/2 results in’ the greater stability of 
the Dufort-Frankel scheme. Even though terms involving S at 
time t + 2At appear on both sides of the equation, implying 
an implicit system, the term on the right may be transposed 
to the left - such schemes are called semi-implicit or 
pseudo-implicit as they may be rewritten in an explicit 
form. In this case, (3.4.5) becomes 

S(J,n+2) {1+2r} = {1-2r} S(J,n) + 2r {S(J-1,n+l) + S(J+1,n+1)}(3.4.6) 
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Numerous difficulties arose when this method was 
applied to the diffusion terms (b) in (2.4.11). As the time 
step was reduced, various parameters of the circulation 
(e.g. the extremum of the streamfunction field) were noted 
to behave as if the truncation error associated with the 
scheme was very large. Note that the Dufort-Frankel scheme 
hes htroncationn errormoft.Od Atte tehy*) +e At /Ayor os. 

To examine this more closely, another formulation was 
applied which is only first order in time. This scheme uses 
a forward-time, centered-space finite-difference technique 
and was applied to the prediction equations. In terms of 


(3.4.4), it can be expressed as 


unas il S(J,n) _ x, son) a eee 27) 


Again, with r = K At/Ay?, this becomes 

S(J,ntl) = r S(J-1,n) + {1-2r} S(J,n) + r S(J+1,n) ‘ (3.4.8) 
The stability criterion for this representation, as noted in 
TaDLe7 3125-795 


KeAt 


af 
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As the time step was reduced using (3.4.8), more consistent 
results were obtained than with (3.4.6), even though the 
latter was second order in time. This revealed the existence 
of a numerical problem when the Dufort-Frankel scheme was 
implemented and incited further exploration. For this 
purpose, a simple one-dimensional heat diffusion (or 


conduction) equation was considered, for which analytical 
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solutions were available. Initial and boundary conditions 
Similar to those of the valley model under consideration 
were applied and numerical results using both the Du- 
fort-Frankel and the forward-time, centered-space schemes 
are presented in Section 3.5 for comparison. Further dis- 
cussion at that time is made regarding the proper choice 


of a diffusion scheme for the prediction equations. 


3.4.4 Other Term (c) 

The solenoid term, denoted by (c) at the beginning of 
this section, expresses the non-linear dependence of & and 
6. The relationship is such that a theoretical investigation 
of the behavior of the various possible numerical schemes is 
precluded. Accordingly, through accuracy considerations 
alone, the centered-time, centered-space formulation was 
selected. The truncation error for this scheme is 
OCAt’ + Ay-). This formulation was usec by Thyer (1966) and 
Orville (1964, 1965, 1968) in their models of flow above 
Sloping terrain and is recommended by Mesinger and Arakawa 
(a7 Cua. 

An example of the finite-difference representation 
applied to 


oe & 28 (3742.10) 
at 6 oy 


is 

E(J,K,n+2) - E(J,K,n)_ g {6(J+1,K,n+1) - 6(J-1,K,n+1)}(3,4.11) 
2At 6(J,K,n+1) 2hy 

Using this scheme, changes in & over the time period from t 


to t +*2At arev evaluated at the middle of the) interval’ and, 
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for this reason, it is often called the mid-point rule. 


3.5 One-dimensional Heat Conduction Equation 

As indicated in Section 3.4.3, analytical and numerical 
solutions to the one-dimensional heat conduction equation 
were obtained in order to examine the accuracy of the numer- 
ical solutions. This is the classical partial differential 
equation of mathematical physics used to describe the 
conduction of heat in a solid body. Here, it is applied to 
predict the potential temperature, 6(z,t), determined by 
diffusion away from a steadily cooling surface located at 
z = 0. The top of the atmosphere is assumed to be at a 
height z = H. The partial differential equation governing 


this transfer process can be written 


96 . 976 
ve oe rradahy Ur aens Vet IA (yn Seal) 


It is assumed that the temperature at the top of the atmos- 
phere remains constant, whereas that at the surface 15 a 
linearly decreasing function of time. Therefore, the boun- 
dary conditions are 

OCH, &)iaeT, t>o0 ; Gay 5e02)) 

6(0,t) = 6(0,0) - at, teon0, (32523) 
where a = constant surface cooling rate. Initially, some 
vertical eee paeion of temperature is assumed so that 


6(z,0) = f(z), 0, <). 2: 4H. 634544) 
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3.5.1 Analytical Solution 
Ananalytical solution to (3.5.1), subject to the 
non-homogeneous boundary conditions, (3.5.2) and (3.5.3), 
and initial conditions 
0(z,0) = f(z) = (T - 6(0,0))z/H + 6(0,0), Oneeze<cH 54 BS. 55) 


is derived in Appendix C to be 


Q(z,t) = } { (-2H20/ (nm) 3K) exp(-n*n2Kt/H?) sin(n7z/H} 
n=1 


Qa 
+ — (z3/6H - 22/2 - zH/3) + (8(0,0) - at) (1 - z/H) (3.5.6) 
K 


+ 2T/H, Or<eZe<s Heute oF 0 


Attempts at solution with 6(z,0) satisfying the Poisson 
equation, (2.2.5), as in the valley model, failed and a 
Simple linear temperature profile was assumed as an 
approximation. The two initial fields differed at most by 
020008 "GC". 

The potential temperature fields given by (3.5.6) are 
plotted for various times in Figures 3.6 and 3.7. The 
abscissa gives the total cooling relative to t = 0 seconds. 
Here, the infinite series in (3.5.6) was truncated at 500 
terms. After 20,000 seconds or = 333 minutes, the tempera- 
ture profile varied almost linearly from the cooled surface 
temperature, 6(0,t) = 6(0,0) - at, to the constant tempera- 


ture at the upper boundary, @(H,t) = T. At earlier times, 
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HEIGHT (MN) 


=017) “ote -0.5 “0.4 -0.3 -0.2 O01 0.0 
TEMPERATURE CHANGE (C°) 


Figure 3.6 Potential temperature deviation from initial 
State versus height for one-dimensional diffusion equation 
(K = 1.0 m? s-') out to 1200 seconds 
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HEIGHT (M) 


-6 ae 
TEMPERATURE CHANGE (C2) 


Figure 3.7 Potential temperature deviation from initial 
State versus height for one-dimensional diffusion equation 
(Kieu =1.0 m* s-') out, to 20000, seconds 
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e.g. Figure 3.6, the non-linear variation with height is 
very evident. Gradually the cooling specified at the surface 


a@iffuses upwards to greater and greater heights. 


3.5.2 Comparison with Numerical Solution 

Generally, the two numerical schemes examined (Du- 
fort-Frankel and forward-time, centered-space) compared 
favourably with the analytical results. A sample computation 


aSuciven in Tabley3.3, upito 20.m height at) t 2000 sec- 


@nds, using Az =.2,.9, Ot = 0.5,seconds,, and) K ee Ouemiae S ous 
These values satisfy the stability criterion for the 
LoOrward-time,. centered-space formulation, (3.4.9). At z= 0} 
the numerical solutions are exactly the analytical ones, as 
specified bythe poundary) condition, (3.5.3). “Poruayconstan: 
grid size, the numerical solutions tended to overshoot the 
analytical ones, as much of the inaccuracy is involved in 
Az, the grid size. However, as both At and Az are reduced, 
both numerical methods approach the true analytical one. 
Thus, it appeared that either numerical scheme would behave 
Satisfactorily when applied to a simple purely diffusional 
problem. 

At this point, the numerical and analytical results 
were compared with the two-dimensional valley model results. 
Initially, temperature changes are brought about, according 
to (22252), virtually as a result icf vertical diffusion from 
the cooling ground surface alone, Since both horizontal dif- 


fusion and advection are small. The valley model predicted 
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Table 3.3 Comparison of analytical solution to one-dimen- 
Sional diffusion equation with numerical solutions obtained 
using Dufort-Frankel and forward-time, centered-space 
schemes after 2000 seconds 


Analytical Solution | Numerical Solution |} Numerical Solution 


0 282 .8738034 282 .8738034 282 .8738034 


283.0376210 283.0376174 283 .0376198 


283.1930005 


283 .1929934 283.1929983 


283.3403911 283 .3403806 283 . 3403878 


283.4802385 283 .4802247 283.4802342 


283.6129835 283.6129667 283 .6129783 
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temperature profile above the valley trough (J = 1) was 
selected for comparison since this column is most remote 
from the sloping surface and would, thus, have the smallest 
contribution from horizontal diffusion. Table 3.4 compares 
the total cooling predicted (relative to t = 0 seconds) out 
as far as 200 seconds from the various numerical diffusion 
schemes applied to the valley model with the analytical 
one-dimensional solution. The valley model results shown are 
for an integration with Ay = 40 m, Az = 4m, and At = 2 sec- 
onds. 

These results are consistent with those previously 
described, i.e. the various numerical diffusion schemes 
perform adequately when used on a purely diffusional problem 
as in the early stages of the evolution of the model when 
advection processes are quite small. The vorticity equation, 
(2.4.11), however, is a different matter. As soon as surface 
cooling commences, the solenoid term (g/@ 06/ay) becomes 
important and the vorticity changes are not solely the 
result of diffusion of vorticity. When the two different nu- 
merical diffusion schemes are applied to (2.4.11), certain 
discrepancies arise. 

in) Tables 3. 5 and=3.26; thes grid point values of fy 
Te Of /ot| es and of/dt| || are given for each of the 
numerical schemes at grid points located in a vertical line 
midway up the slope. These are for a valley model integra- 
tion out to 20 seconds, with input parameters of Ay = 40 m, 


Az.=>4.m,, and. Ai =.4) seconds... Theypattern of. vorticity 
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Table 3.4 Comparison of amounts of cooling after 200 seconds 
obtained from analytical solution to the one-dimensional 
diffusion equation and valley model solutions obtained using 
Dufort-Frankel and forward-time, centered-space schemes 


Height | Analytical Solution ee Model | Valley Model 
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Table 3.5 Vorticity changes predicted by various terms in 
the valley model out to 20 seconds using the Dufort-Frankel 
scheme 
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Table 3.6 Vorticity changes predicted by various terms in 
the valley model out to 20 seconds using the forward-time, 
centered-space scheme 
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changes due to diffusion shows quite a difference between 
the two schemes. The solenoid terms, which depend only on 
the temperature fields, are very Similar as the two tempera- 
ture evolutions are not drastically different. However, 
after only 20 seconds, the vorticity resulting from the Du- 
fort-Frankel formulation is more than 28 percent higher than 
that resulting from the forward-time, centered-space scheme, 
and 40 percent higher after only 52 seconds. This 
discrepancy can be traced mainly to the different techniques 
involved in diffusing the vorticity values. 

The Dufort-Frankel scheme is a three-level scheme in 
time, and, as a result, the changes at time t depend on 
values at t - At and t - 2At. If these values have 
previously been altered only through diffusive processes, 
then reasonable results are obtained. With the solenoid term 
coming into play at every time step, however, the values at 
the intermediate time,;”t'— At, “cause a’shock or: disturbance 
in the diffusion computation. The forward-time, 
centered-space scheme, on the other hand, considers only the 
field at a single previous time step and results in 
reasonable diffusion changes over that time interval. 

Since the problem is suggested to lie in the magnitude 
of the time step, the scheme which 1S most nearly correct is 
expected to show the least difference when the time step is 
reduced. This turns out to favour the use of the 
forward-time, centered-space formulation. Table 3.7 shows 


the vorticity values predicted above the mid-slope region 
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Table 3.7 Vorticity values (s~') predicted above mid-slope 
region after 200 seconds using the Dufort-Frankel and 
forward-time, centered-space schemes and various time steps 
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after 200 seconds using the two numerical schemes and 
various time steps. Clearly, as the time step is reduced, 
the results of the Dufort-Frankel scheme seem to approach 
the results of the forward-time, centered-space scheme, as 
anticipated. This is also reflected in Figure 4.1, which 
depicts the evolution of the minimum streamfunction value 
with time for the two schemes. Based on these results, the 
latter finite-difference formulation was chosen to 
approximate the diffusion terms in this model. This scheme 


was used by Pielke (1974). 


3.6 Liebmann Sequential Relaxation Technique 

A solution to the equation relating the streamfunction 
Vearueseevicr hk, tos tO tne VOrticity! values Ee (gek, t)= 
(2.5.3), was obtained numerically using Liebmann sequential 
relaxation’ (Haltiner@and Wi bitamsa (71980) aiProme (2.4.11), 
ehemvorcicity field’ was) predicted) forvallsgrid points. At 
each time step, the velocity components were required by the 
model (for use in predicting contributions to 00/dt and 
d&/dt through advection) and these were obtainable from 
(2.5.1) and (2.5.2), once the streamfunction field 
corresponding to the vorticity field was known. Since the 
streamfunction value everywhere on the boundary was set to 
zero, the relaxation technique was applied only on the 
interior of the grid away from the boundary. Details of the 


technique follow. 
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Let the nth guess at the streamfunction value at grid 


point (J,K) be y°(d,k). Then, 


E(J,K) + V? y"(J,K) = R'(J,R) (3.6.1) 
where R"(J,K), called the residual, is a measure of the 
error of the guess. If w"(J,K) satisfies (2.5.3) exactly, 
then R"(J,K) will be zero. Since a guess at w(J,K) is 
unlikely to be the correct solution, the residual, R"(J,K), 
will not be zero at every point of the grid, if anywhere at 
all. However, successive estimates of the streamfunction 
field are made, according to the algorithm below, until all 
residuals are less than some preassigned tolerance level, e. 

Expanding (3.6.1) into its centered finite-difference 


representation yields 


E(J,K) + {y" (J+1, K)=2y" (J, K)+y" I-1,K)} 


Ay 
ra PGnZive” Le K)+ "GK-1)} = RG K) 
Az 


We wish to determine a new guess, wW"*'(J,K), which will 


Cor Gee) 


reduce the residual to zero. Therefore, 


EC) Ket (y® (341K) 29", K-40" (J-1,K)) 


by 
+1 
+ {py (J, K+1)-2y" UK) +9" (J, K-1)} | 0 
Wy? 


Subtracting 5.6.5). from (3.6.2). eancusolvingatorethesnew 


(Cenc) 


Guess ewie. ( Jo Kk) results. in 


yrtt sx) = y"(s,K) + R7C,K) 1 
iGe Sa) (3°64) 


This shows that the (n+1)th guess of w(J,K), which is 


increased over the previous guess by the quantity 
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R°(J,K)/(2/Ay? + 2/Az*), will reduce the residual to zero at 
that particular grid point. However, such a correction at a 
particular point affects the residuals at adjacent grid 
points. Nevertheless, each guess made according to (3.6.4) 
represents overall progress in reducing all the residuals, 
and the method is guaranteed to converge towards the true 
streamfunction values. 

It has been noted that an overrelaxation technique will 
result in faster convergence. In this case, a larger 
correction than that given by (3.6.4) is added to the old 


guess, using 


n+l n n 
y (J,K) = p (J,K) + a RGi,K) ’ 
(Tay2 + 2/422) URES (3.6.5) 


where a 1S an overrelaxation factor and usually 1 <a < 2. 
Theoretical estimates of aw have been made based on the grid 
size but numerical experiments are sometimes used to 
determine the optimal overrelaxation factor which produces 
fastest convergence. 

The technique of sequential or Liebmann relaxation was 
usually found to converge even more rapidly. With this 
method, each new guess made using (3.6.5) is immediately 
incorporated into the determination of the next new guess at 
adjacent grid points. As well, this scheme is more 
economical with respect to computer Storage. It is this 
scheme, with a = 1.2 (determined experimentally), which was 
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3.7 Summary of Method of Solution 

A summary flowchart which describes the numerical 
schemes utilized in the valley model is given in Figure 3.8. 
Certain peculiarities not previously mentioned should be 
noted. For the first increment in time, forward or Euler 
differencing must be used as the values at the midpoint of 
the time interval are unknown. Subsequently, at least for 
the solenoid term, a three-level scheme is possible. 

According to the stability criterion for the 
forward-time, upstream-space formulation of the advective 


terms, (3.4.2), the time step used must satisfy 
Atts a anda At <Pae Ces IRA) 


As v and w change with time, an estimate of their maximum 
values is used in (3.7.1) to determine the corresponding 
maximum time step. As the magnitude of the circulation 
increases, provisions are made in the computer code to 
ensure that (3.7.1) remains satisfied. An insurance factor 
of 0.80 was used so that the time step required 
oversatisfied the above inequality. If the time step became 


no longer suitable, it was halved successively until both 


at and at < SE Ae G8 Weer 2) 


were true. The first increment in time after such a 
reduction would again require forward differencing as in the 
beginning. This procedure is illustrated in Figure 3.9 for a 


case with an initial time step of 5 seconds. It is 
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Figure 3.8 Flowchart of finite-difference schemes used in 
valley model 
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Figure 3.9 Example of time differencing used in the valley 
model (F = forward, C = centered) 
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postulated that at t = 15 seconds, the time step must be 
reduced to 2.5 seconds in order to satisfy (3.7.2) and the 
appropriate type of differencing (either forward or 


centered), 21s indicated. 
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Chapter 4 


RESULTS 


4.1 Introduction 

In this chapter, the results of numerous valley model 
integrations are presented. Initially, the input parameter 
values for a standard run are considered. Sensitivity 
analyses are presented in which several of these parameters 
are varied, one at a time, in order to determine the effect 
of each. Some discussion, of a numerical nature, iS given in 
regard to the choice of time step suitable for a particular 
run, etc. along with the actual CPU costs. 

The results are presented in a variety of manners : 
two-dimensional streamline plots, time and space evolution 
of the streamline center, development of slope winds with 
time, cooling amounts at different heights in different 
parts of the valley, relative contributions of advection and 
diffusion, ate Variations in the geometry of the grid 
configuration (location of upper boundary, presence of hori- 
zontal valley floor versus 'V-shaped region, angle of slope) 
are considered. Sensitivity to the vertical diffusion 
coefficient under an initially stable regime is provided 
through two integrations with quite different values. 
Several initial stabilities are implemented, with vertical 
temperature gradients ranging from -0.04 C°/m up to 


+0,.01°¢2%/mu( the yvnevtral vease): 
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Finally, the results from simulations in which values 
most appropriate to the microclimatology of the river valley 
are presented. Limited comparison of the model with 


observations made within the valley is carried out. 


4.2 The Standard Run - Run 1 


4.2.1 General 

In this section, the results from what was considered a 
Standard run are presented. As Run 1 will be frequently used 
for comparison with other model integrations, analysis of 
the results obtained is presented in greater detail than 
for some subsequent runs. The input parameters used for 
the standard run are given in Table 4.1. Those variables 
denoted by an asterisk were allowed to vary as part of the 
sensitivity analyses and these results are shown in later 
Sections. An initial valley trough temperature, T(1,1,0), of 
2/0. Zoeele WSO SeUlS COs Neat ASSO ueOlaGUNSoncLOndawithwant initial 
valley trough pressure, P(1,1,0) of 930 mb. Using these 
values and the various lapse rates, the initial potential 
temperature profile was computed according to Section 3.3.1. 
Under stable conditions, isotherms were assumed to be hori- 
zontal within the valley cross-section, as shown later in 
Figure 4.7. Table 4.2 provides a summary of parameter values 
used in subsequent runs. 

Preliminary tests were conducted in order to determine 
the appropriate time step. AS an estimate of what the time 


step should be, the stability criteria for the various 
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Table 4.1 Parameter values for valley model in the standard 
run (Run 1) 


CK), 7 &), 
(Kea) oem (Keds a ® 
vertical ei size (below 72 m) 4m 
(above 72 m) 12 m 
PHYSICAL GRID PROPERTIESfhorizontal grid size 40 m 


height of upper boundary * 108 m 


inclination of slope * .09978=5.7° 


time step * 
INITIAL CONDITIONS 
lapse rate 


fe) 


* allowed to vary for sensitivity analysis 


Table 4.2 Summary of parameter values for Runs 2 - 8 


* except upper boundary at 144 m 


* except V-shaped valley cross-section 


* except angle of inclination of slope (a) = 0.2 


* = = 2 Syl 
except Ke (KK), 0.01 m* s 


* except initial lapse rate = ane = +0.01 C°/m 


-0.04 C°/m 


* except initial lapse rate 


* except initial lapse rate 


B/C» surface cooling 


rate = 2.5 C°/hr, (K), = (Kp), = 0-0 


* as in Run 1 (see Table 4.1) 
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numerical schemes were examined. For the advective terms, a 
time step of about 20 seconds or fess was needed, based on 
(3.7.2) and considerations outlined in Section 3.7. The dif- 
fusion scheme chosen was more restrictive, however, at least 
in terms of the vertical grid spacing. Substituting into 
(3.4.9), and recalling that integration was in a forward 
sense over two time steps, revealed the requirement of a 
time step of 4 seconds or less. 

The location of the minimum streamfunction value at a 
given time showed little variation as the time step was 
reduced. For this reason, the magnitude of the streamfunc- 
tion minimum, typically located somewhere above the sloping 
ground, could be compared under the various numerical 
schemes and time steps implemented. In Figure 4.1, the mini- 
mum streamfunction value predicted is plotted versus time 
for various time steps. The problematic Dufort-Frankel 
scheme is also presented for comparison. 

The forward-time, centered-space diffusion scheme 
appears to converge consistently to a single value, and, at 
first glance, a time step of 4 seconds would seem 
sufficient. However, upon closer examination, an oscillation 
in sign appeared in the diffusion changes predicted at ver- 
tical grid points above the valley trough when a time step 
of 8 or 4 seconds was used. An example of this instability 
US SHOWN i nshigune 4.2... In which, thevamount..of .cooling 
predicted after various times 1s plotted versus height above 


the valley trough usSing a time step of 4 seconds. A similar 
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Figure 4.1 Minimum streamfunction versus time for Run 1 
using various diffusion schemes and several time steps 
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Figure 4.2 Cooling predicted by valley model versus height 
uSing forward-time, centered-space diffusion scheme and a 4 
second time step 
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oscillation occurred at an earlier time when a time step of 
8 seconds was used, again only above the valley trough. 
Above the sloping surface, it did not appear and is 
therefore not evident in Figure 4.1. With a time step of 2 
seconds, such an oscillation did not: occur, at least out to 
integration times where other difficulties arose and it was 
chosen for this reason. 

A compromise between the cost of each integration using 
a certain time step and the results achieved was made. This 
run, with a time step of 2 seconds, required 1000 seconds of 
CPU time in order to integrate out to 2000 seconds. A 
Similar run, with a time step of 4 seconds, required 408 
seconds of CPU time to reach 1500 seconds, indicating the 
non-linear relationship between time step and CPU 
requirements. Most of the CPU time is used up in the 
relaxation process. With smaller time steps,smaller changes 
in the vorticity are produced and thus fewer relaxations are 


required in order to achieve a certain accuracy. 


4.2.2 Circulation Development 

Initially, negative vorticity was generated close to 
the ground surface through the solenoid term, according to 
(2.4.9). This resulted in a corresponding downslope flow 
close to the sloping surface with a weaker return flow 
aloft. Figure 4.3 depicts the streamlines predicted after 
600, 1200, and 1400 seconds. Minimum streamfunction values 


and their location are noted. Out to an integration time of 
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Figure 4.3 Streamfunction contours (Aw = 2.0 m? s-') att 
600, 1200, and 1400 seconds for the standard run (Run 1) 
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about 1250 seconds, the qualitative nature of the flow 
pattern remained essentially the same. There was a 
well-developed downslope wind, which was confined to a thin 
layer near the surface. A weak updraft existed above the 
valley trough which spread horizontally at some height. In 
the upper parts of the region, there exists this mainly hor- 
izontal flow away from the trough, with weak downward motion 
towards the ridge. 

Figure 4.4 indicates the trajectory of the streamline 
center obtained a.)from the available grid point values, and 
b.) from a non-linear interpolation procedure. This method, 
although somewhat arbitrary, was felt to offer a significant 
improvement in the estimation of the location of the stream- 
line center. From this Figure, the streamline center was 
generally seen to move up the sloping surface at some height 
above it (24 - 28 m). It originated near the middle of the 
Slope, strengthened with time, and eventually moved up and 
Slightly further away from the surface. Little movement of 
the center occurred once the ridge height was reached. 

Along the sloping ground surface, the velocities 
predicted by the model were quite uniform. The rate of 
development of the downslope wind appeared to decrease 
beyond about 1000 seconds of integration time, as it 
approached a speed of about 50 cm s~'. An illustration of 
the relationship between the streamfunction field and the 
velocity components themselves iS given in Figure 4.5. 


Figure 4.5a consists of the isopleth contours of the 
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Figure 4.4 Trajectory of streamline center for the standard 
run (Run 1) obtained from (a.) available grid point values 
and (b.) interpolation. Numeral indicates time in hundreds 
of seconds 
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Figure 4.5a Circulation predicted after 400 seconds for the 
Standard run (Run 1) - streamfunction (Ay = 0.5 m? s-', 
dotted) and velocity (A|V| = 0.01 m s-', solid) contours 
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Figure 4.5b Circulation predicted after 400 seconds for the 
Standard run (Run 1) - horizontal velocity component 
Uive=" 0.04 m,s-!) contours 
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Figure 4.5c Circulation predicted after 400 seconds for the 
Standard run (Run 1) - vertical velocity component 
(Aw = 0.005 m s~') contours 
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streamfunction after 400 seconds, along with contours of the 
speed of the corresponding flow given by ha = ery sews he? 
Large velocities were present between the surface and the 
Streamline center where the gradients in streamfunction were 
high. Surface velocities at this time reached 11 cm s~'. 
Figures 4.5b and 4.5c illustrate the spatial pattern of the 
two velocity components, v ana w, respectively. Figure 4.5b 
is quite Similar to Figure 4.5a as the vertical velocity 
components were relatively small. Since the gradients in 
temperature were large under this stable regime, however, 
these small vertical velocities were able to alter the tem- 
perature fields significantly through vertical advection. 
Initially the maximum upward velocity occurred not along the 
Symmetry axils but rather several grid points to the 
interior, along J = 3, above where the sloping surface 
intersected the horizontal valley bottom. Gradually, the 
maximum vertical velocity moved towards the valley trough 


boundary. 


4.2.3 Evolution of Thermodynamic Fields 

The thermodynamic field is altered through both ad- 
vective and diffusive/radiative processes. The winds 
generated act to bring warmer air down from above and vice 
versa. Diffusion acts to mix the cooled air near the surface 
with the warmer air aloft. The relative magnitudes of these 
processes change with time and with position in the valley 


and are given in Table 4.3 for a grid point 4 m above the 
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middle of the slope (at J = 8, K = 7) at 200 second 
intervals. Each term in (2.2.2) was evaluated using the 
various numerical schemes outlined in the previous chapter. 

Horizontal advection was found to always act as a 
cooling process. The horizontal wind component was always 
negative or directed away from the cooling slope, and 
therefore acted to transport cooler air away from the 
surface. The vertical component was also negative and 
approximately equal to 10 percent of the horizontal 
component at this grid location. The vertical temperature 
gradient was large and positive here, resulting in the 
transport of warmer air from above towards the surface. At 
all times, the amount of vertical advection exceeded that of 
horizontal - a layer of warming through advection was the 
overall result just above the surface. 

Referring again to Table 4.3, both horizontal and ver- 
EPrcaly Q’LfLusion terms resulted in scooling at. this guid 
location. Due to the large cooling rates at the surface, the 
temperature gradient between the surface and this grid point 
exceeded that between this grid point and the next one away 
from the surface. Thus, the Laplacian in temperature 
(d/dy 06/dy + 3/dz 06/dz) was clearly negative. Since the 
horizontal and vertical coefficients of diffusivity were 
equal and Az was much smaller than Ay, cooling through ver- 
ticai diffusion greatly exceeded horizontal diffusion. 

Figure 4.6 illustrates the relative contributions of 


advection and diffusion in various parts of the valley after 
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Table 4.3 Contributions to 06/dt from advection and 
diffusion just above mid-slope (J=8,K=7) for Run 1 at 200 
second intervals 
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400 seconds (Figure 4.6a) and after 1200 seconds (Figure 
4.6b). Generally, warming through subsidence occurred along 
the slope and above the valley ridge. Greater cooling 
through diffusion occurred at these locations, however, but 
as the velocities increased with time the difference in the 
magnitudes of advection and diffusion diminished. Cooling 
through diffusion was a maximum just above the ridge surface 
where cooling below and warming above together created large 
differences in the vertical temperature gradient. 

The actual potential temperature fields predicted can 
be examined in a number of ways. A sequence in time of 
two-dimensional plots of potential temperature over the val- 
ley cross-section are shown in Figure 4.7 at 400 second 
intervals. Early in the model integration, distortion of the 
horizontal isotherms through cooling at the surface was very 
evident. Above the ridge, slight packing of the isotherms 
appeared as a result of cooling from below and warming 
through subsidence from above. This acted to strengthen the 
vertical temperature gradient close to the ground. The ver- 
tical temperature gradient was 0.014 C° m™' at t = 0 sec- 
onds. After 1200 seconds, the potential temperature at the 
ridge surface had dropped to 284.046 °K, whereas that 32 m 
above had increased to 285.170 °K - a vertical temperature 
gradient of 0.035 C° m-' resulted, over twice the initial 
value. At the trough axis, there was cooling at the surface 
and somewhat less cooling aloft - a vertical temperature 


gradient of 0.022 C° m-'resulted after 1200 seconds. The 
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Figure 4.7a Potential temperature field at t = 0 seconds for 
the standard run (Run 1) 
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Figure 4.7b Potential temperature field at t = 400 seconds 
for the standard run (Run 1) 
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Figure 4.7c Potential temperature field at t = 
for the standard run (Run 1) 
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Figure 4.7d Potential temperature field at t = 


for; the. standara run. (Run: 1) 
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Stability had increased more at the ridge height than at the 
trough level. 

Another useful manner in which to examine the tempera- 
ture changes as functions of time and position is through 
plots of the vertical temperature profile in different parts 
of the valley cross-section for several times. Figures 4.8, 
4.9, and 4.10 illustrate vertical profiles of the cooling 
predicted as a function of height at the trough, mid-slope 
and ridge positions, respectively. Comparing Figure 
428 with Figure 3.6, in which the amount of tooling 
predicted through vertical diffusion alone 1s given, reveals 
the contributions which advection and horizontal diffusion 
make to the temperature deviation. Figure 4.11 provides such 
a comparison. The difference between the valley model 
results and the analytically derived one-dimensional . 
solution are given at several times in this set of curves. 
At early times (up to 200 - 300 seconds of integration 
time), there was very little difference between the two 
solutions, indicating that the predicted valley model 
changes occurred through vertical diffusion, as expected. 
The difference between the two solutions was zero at the 
upper and lower boundaries, where exactly the same boundary 
conditions were specified in both cases. Close to the ground 
surface, vertical diffusion played a larger role in altering 
the temperatures than in the atmosphere above, and so the 
two solutions differed less and less as the surface was 


approached. 
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Figure 4.8 Potential temperature deviation from initial 
state versus height above trough boundary (J=1) for the 
Standard run (Run 1) 
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Figure 4.9 Potential temperature deviation from initial 
state versus height above mid-slope region (J=8) for the 
standard run (Run 1) 
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Figure 4.10 Potential temperature deviation from initial 
State versus height above ridge boundary (J=16) for the 
standard run (Run 1) 
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Figure 4.11 Difference in cooling predicted by valley model 
in Run 1 and by one-dimensional vertical diffusion equation 
(Kumui.0 m- S-'))as au functions ofeneiagnt 
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Beyond about 1260 seconds, the temperature profile 
above the trough, shown in Figure 4.8, began to exhibit some 
unusual features. The source of this extreme warming aloft 
was revealed upon re-examination of the circulation regimes 
depicted in Figure 4.3 at this time. A counterclockwise cell 
had begun to develop well above the valley floor just after 
1200 seconds. This resulted in downward motion along the 
trough boundary, with accompanying warming through subsi- 
dence. The origin of this counterclockwise motion was the 
larger vertical velocities which occurred along J = 3 
compared to those along J = 1, as is evident in Figure 4.5c. 
A maximum upward velocity of about 0.9 cm s~' existed along 
J = 3 at 400 seconds whereas the maximum along J = 1 (the 
lateral boundary) was about 0.8 cm s~' at that time. After 
1200 seconds, these velocities had increased to 3.5 cm s~“' 


anand<;Stemns”! 


along J = 3 and J = 1, respectively. The 
resulting differential in the amounts of vertical advection 
created a slightly cooler area at some distance from the 
boundary. As the magnitude of the vertical velocity 
increased, the amplitude of the wave in the temperature 
field close to the trough increased. Through these 
variations in the horizontal temperature gradient, the vor- 
ticity eventually gained positive values well above the val- 
ley trough and this corresponded to a counterclockwise flow 
in this area beyond 1260 seconds. 


Figures 4.9 and 4.10 illustrate the vertical tempera- 


ture profiles above the mid-slope region (J = 8) and ridge 
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(J = 16), respectively. The general shape of the curves in 
Figure 4.9 is similar to those in Figure 4.8, especially at 
early times. Later, warm advection along the slope acted to 
retard the effect of the cooling ground surface below. In 
Figure 4.10, the increase in potential temperature at 
heights of 30 m or more above the ground is clearly shown. 
This was due to, aS mentioned previously, the downward flow 
above the ridge which created warming through subsidence. At 
low levels, cooling persists due to mixing of the air near 


the cooled surface. 


4.3 Effect of Location of Upper Boundary - Run 2 

A second model integration was performed using the same 
Pipiieoalame vel osaswinenun le ( listed sine labled i excepr 
that the upper boundary was located at a height of 144 m 
rather than 108 m. Thus, a vertical grid spacing of 4 m was 
used from the surface up to 72 m, and a Spacing of 12 m from 
72 m up to 144 m. An initial time step of 2 seconds was 
again used, with integration out to 1400 seconds requiring 
706 seconds of CPU time. 

The predicted streamfunction fields were found to 
exhibit very similar properties to those of Run 1, at least 
qualitatively. The return flow aloft reached to somewhat 
greater heights but because of weaker vertical gradients 
between the streamfunction center and the upper boundary, 
smaller horizontal winds resulted. For example, at 600 sec- 


onds, the maximum horizontal wind was about 5.1mm s"' at a 
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= when 


height of 72 m in Run 1 but was only about 3.6 mm s 
the flow was permitted to greater heights. Vertical winds 
near the trough boundary were slightly raised over those of 
Run 1, with the result that the counterclockwise cell 
appeared somewhat earlier following the reasoning given in 
Section 4.2.3. Development of downslope winds near the 
sloping surface was generally identical. The trajectory of 
the streamline center appeared to move farther up and away 
from the slope than in Run 1, eSpecially within the early 
Stages of the integration. Significant movement of the 
center ceased once the height of the ridge was reached, as 
2nisRun 16 

The temperature profiles out to 1000 seconds or so were 
virtually identical to those for Run 1 for the valley region 
below ridge height. Further above, the greater height of the 
upper boundary permitted temperatures there to be altered 
through both advection and diffusion, although the warming 
Or cooling predicted was quite small. Beyond about 1000 sec- 
onds, large differences existed between the two runs, 
principally due to the cellular nature of the flow in the 


second run by this time. 


4.4 Effect of Absence of Horizontal Valley Floor - Run 3 

In order to determine the influence of including a hor- 
izontal valley floor, one integration within a ‘V-shaped’ 
valley cross-section was performed. All input parameters 


were the same aS in Run 1, except that the two columns of 
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grid points above the flat valley bottom were eliminated. 
This effectively moved the lateral trough boundary over to 
where the sloping ground intersected the previously included 
valley bottom. Again, a 2 second time step was specified 
initially, and 506 seconds of CPU time were used in reaching 
1600 seconds of model time. 

In the earlier stages of this run, out to about 1000 
seconds of integration time, the primary difference was 
found close to the trough boundary, as expected. Maximum 
vertical velocities above where the sloping ground met the 
horizontal are given in Table 4.4 for Run 1 (J = 3) and Run 
3 (J = 1). This enables the same relative positions with 
respect to the sloping ground to be compared. The height at 
which the maximum vertical velocity was found, based on the 
available grid points, 1s also provided in this Table. In 
Run 3, the flow is forced to rise very close to the slope, 
resulting in larger horizontal gradients in the streamfunc- 
tion (and therefore larger vertical velocities) in this 
area. 

The potential temperature profile above the valley 
trough reflected the larger vertical velocities as greater 
cooling through adiabatic expansion was evident. At a height 
of. 52..m, for “examplepethe=resultseforeRunes-indtcated a-pot- 
ential temperature of 284.64 °K after 600 seconds whereas 
Run 1 predicted 284.66 °K. After 1000 seconds, Run 3 and Run 
1 predicted temperatures of 284.46 °K and 284.50 °K, res- 


pectively, at 52 m. The mid-slope temperature profiles were 
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Table 4.4 Maximum vertical velocities above valley trough 
Porerun i (J=3)° and. for Run sy (J=1)).. Quantityesin 
parentheses denotes height at which maximum occurred 
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quite Similar, however, with slightly less warming at low 
levels, and slightly more cooling at upper levels through 
advective differences. Above the ridge, little difference, 


if any, was evident in the temperature profiles. 


4.5 Effect of Steeper Slope - Run 4 

In order to determine the effect on the evolution of 
the circulation of a lower boundary which was twice as steep 
as that of the standard run, a fourth integration was 
performed. Apart from interest in discovering how the model 
would behave under such conditions, this run was felt to be 
useful due to the asymmetric topography of the North 
Saskatchewan River valley. From Figure 3.1, it is evident 
that the valley is much steeper on the north bank as the 
river meanders eaStward cutting into the bank on the outside 
Gorner.. These, arsVope anglevot about ten '(0.22) orp 0.23 
radians would be more representative. For these reasons, the 
slope angle for Run 4 was set to tan-'(0.20) or about 0.20 
radians. As the permissible grid Spacings are closely 
related to the slope angle through (3.2.1), the horizontal 
Grids size) Ay; wWwas*reduced to 20 mo for this: run. Both Az, 
and Azz remained 4 m and 12 m_, respectively. The region 
being modelled here is shaded in Figure 4.12, and is shown 
relative to the region under consideration in the standard 
run. A time step of 2 seconds sufficed since the vertical 
grid spacing, which was the limiting factor in the standard 


run, remained unchanged. In order to reach 1400 seconds, 772 
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Figure 4.12 Schematic illustration of modelled region in Run 
4 (hatched), along with that of Run 1 
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seconds of CPU time were required. 

The intensity of the circulation as development 
proceeded can be obtained from Figure 4.13, in which the 
variation of the streamfunction extremum and the maximum 
downslope wind at the surface with time are given for both 
the standard run and Run 4. The steeper slope resulted in 
larger velocities close to the ground, and after about 800 
seconds, downslope winds of about 0.40 m s-' were predicted. 
Thereafter, a quaSi-steady-state regime existed in the cir- 
culation while cooling proceeded at the ground. The term 
qQuasi-steady-state is used here to indicate that the rate of 
change of vorticity, d&/dt, was close to zero at all grid 
points. This means that the three processes which can act to 
alter the vorticity (advection, diffusion, and the solenoid 
term) were in balance. For example, close to 800 seconds, 
the solenoid term acted to generate negative vorticity just 
above the mid-slope region. As well, the vorticity gradients 
were such that transport by the downslope wind acted to 
decrease the vorticity. However, changes due to diffusion 
were poSitive and approximately balanced the two processes 
above. In the standard run, however, a quasi-steady-state 
Period vaid Not Cccurmpricnm to the: formation off aycounter- 
clockwise cell above the trough as in Run 4. 

Increased vertical velocities relative to the standard 
run resulted in increased cooling above the valley trough. 
Similarly, above the ridge, increased warming through subsi- 


dence was predicted for Run 4 as the downward velocities 
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Figure 4.13a Minimum streamfunction versus time for the 
Standard run (Run 1) and for Run 4 
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Figure 4.13b Maximum surface downslope wind speed versus time 
for the standard run (Run 1) and for Run 4 
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were expected to be 2 - 3 times larger than in Run 1. Since 
the downslope winds at a given time were only slightly 
larger in Run 4 than in Run 1, the mid-slope temperature 


profiles were quite similar. 


4.6 Effect of Smaller Vertical Diffusion Coefficient - Run 5 

In order to examine the circulation development when a 
much smaller vertical transfer rate through diffusion was 
used, a fifth integration was executed. With 
Kaos A0f-em*tse pandaky =oi:.0 motsoe,ttameSsteps up £02400 
seconds would, in theory, remain stable, according to 
Criteria based on the diffusion terms (3.4.9). However, the 
horizontal advection term, (3.4.2), required a time step of 
20 seconds or less. For this integration, a time step of 10 
seconds was implemented. This meant that only 247 seconds of 
CPU time were needed in order to reach 2000 seconds of model 
time. 

Through the much reduced vertical diffusion transfer, 
relative to the standard run, the streamline center was kept 
very close to the ground (4 - 8 m away). As well, the rate 
at which the extreme streamfunction value increased with 
time was much less than for Run 1. For example, after 1000 
seconds, the streamfunction attained a value of -2.63 m* s~' 
and was found about 8 m from the slope in Run 5 whereas 
Figures 4.1 and 4.4 indicate an extreme streamfunction value 
of 5.73 -m+ +s" }.to,exisp about i29imefirom the tslopéoin.Run 1 


at that time. However, large downslope winds were the result 
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at the surface due to the large vertical streamfunction gra- 
dient formed close to the ground. Just above the ground 
itself, lower wind speeds than in Run 1 were found. 

Table 4.5 indicates a quantitative comparison between 
the relative magnitudes of advection and diffusion just 
above the mid-slope region. Referring again to Table 4.3, 
which 1S appropriate to Run 1, provides an interesting 
comparison. Horizontal advection acted to cool the air above 
the surface, as in Run 1, but was larger in magnitude. The 
cool air at the surface could not be transferred up and away 
as quickly as in Run 1 because the vertical diffusion rates 
specified were much smaller. Vertical advection was positive 
(warm) and initially smaller than in Run 1 because the ver- 
tical velocities were smaller. After about 800 seconds, the 
vertical gradient in temperature was quite big, resulting in 
large positive vertical advective changes. The sum of these 
two components yielded an overall negative change in potent- 
ial temperature, or cold advection, since the horizontal 
component dominated. This 1s in contrast to the standard run 
in which a layer of warm advection was produced near the 
surface which was able to counteract the cooling through 
diffusion. Horizontal diffusion was small and negative at 
this grid location, but was typically an order of magnitude 
larger than for the standard run. This was due to the slower 
transfer processes in the vertical which created large gra- 
dients in the temperature gradient close to the slope. 


Vertical diffusion was negative and approximately an order 
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Table 4.5 Contributions to 06/dt from advection and 
diffusion just above mid-slope (J=8,K=7) for Run 5 at 200 
second intervals 
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of magnitude smaller in Run 5 than in the standard run. 
Overall, smaller cooling amounts occurred just above the 
ground than in Run 1 even though advection was a generally 


cooling influence. 


4.7 Effect of Stability 


4.7.1 General 

Fleagle (1950) developed a very simple model to 
describe the motion of a compressible fluid due to cooling 
at the bottom by contact with a radiating surface of uniform 
Slope and large extent. He showed that as the air 
accelerated down the slope, adiabatic heating resulted in a 
reversal in the pressure gradient, which retarded the flow. 
As the air decelerated, friction decreased, radiational 
cooling increased the pressure gradient and the cycle was 
repeated. The exact nature of the analytical solution 
depended on the form of the frictional term. If it was 
assumed proportional to the speed, the mean velocity in the 
layer close to the ground varied periodically at first and 
Gradually became constant. The value of this unchanging 
downslope wind was found to be a.) proportional to the net 
outgoing radiation, b.) inversely proportional to the 
thickness of the layer to which cooling extends, and c.) 
inversely proportional to the slope of the ground. 

Tyson (1968) described observations made in 
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downvalley flows up to 6.5 m s~' occurred soon after sunset. 
The wind then steadily weakened throughout the night as cool 
air of marked stability filled the entire valley up to the 
height of the ridge. This showed the variation in the 
Strength of the flow under different stabilities. 

Petkovsek and Hocevar (1971), in an analytical study, 
also included an ambient stratification and were able to 
demonstrate the importance of stratification in determining 
the strength of drainage flows. They found that the more the 
lapse rate in the atmosphere outside the cooled layer 
Giffered from isothermal and the closer it was to adiabatic, 
the larger was the downslope velocity. 

McNider (1982) made a simple extension of Fleagle's 
work in order to examine the impact of various 
stratifications on the oscillatory nature of the resultant 
flow. These results showed that the period of the 
oscillation, as well as the strength of the flow, increased 
with decreasing stability. The periods he found varied from 
20 - 90 minutes as the lapse rate decreased from 8 C° km-@' 
to 2 C° km~'. These periods are somewhat larger than the 
model integration times of the valley model presented in 
this thesis and are not considered here. The mean downslope 
flow varied from 0.06 ms~' to 0.29 m s~-' under the same 


variations in stability, above a slope of about 20 °. 
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4.7.2 Neutral Case Initially - Run 6 

In order to simulate neutral stability conditions ini- 
tially, a model integration was performed using a dry 
adiabatic lapse rate. AS a result, the potential temperature 
field was constant everywhere initially. Accordingly, larger 
downslope velocities were expected since the effect of 
decreasing stratification is to accelerate the drainage 
flow. The stratification within a valley can range from near 
adiabatic at sunset, when drainage flows begin, to many 
degrees per kilometer as cool air pools within the valley. 
This neutral run required only 663 seconds of CPU time to 
reach 2000 seconds of model time. 

In Run 6, uni-cellular flow was maintained out to 2000 
seconds, but some indication was evident of the usual coun- 
terclockwise cell above the valley floor. The rate of 
intensification of the streamfunction center with time is 
plotted in Figure 4.14, along with the maximum downslope 
wind predicted. After about 1600 seconds, the extreme 
Streamfunction value tended to decrease slightly. Downslope 


-' were obtained during this period. Up to 


speeds of 0.9 ms 
about 800 seconds, the predicted maximum slope wind in this 
neutral case did not differ greatly from those of the stan- 
dard run. By 1000 seconds, however, the winds for Run 6 were 
about 25 percent greater than for Run 1. The trajectory of 
the streamline center showed some tendency to move away from 


the slope at an earlier time than in the standard run. From 


1000 to 2000 seconds, little movement of the center 
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Figure 4.14a Minimum streamfunction versus time for the 
Standard run (Run 1) and for Runs 6 and 7 
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Figure 4.14b Maximum surface downslope wind speed versus time 
for the standard run (Run 1) and for Runs 6 and 7 : 
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occurred. Since, from Figure 4.14, the magnitude of the 
Streamfunction center continued to increase, at least out to 
about 1600 seconds, increases in the vertical streamfunction 
gradient (and, therefore, the horizontal velocity) were 
substantial. 

Table 4.6 gives the relative magnitudes of the various 
processes acting to alter the temperatures, valid at a grid 
point which is located 4 m above the slope (compare with 
Table 4.3 for Run 1). Horizontal advection of temperature 
was negative and approximately equalled that for Run 1 up to 
700 seconds or so as the gradients and velocities involved 
were about the same just above the surface. Later, stronger 
velocities produced more cooling by horizontal advection, 
relative to Run 1. Vertical temperature advection was 
positive at this grid location in both cases, but was 
Smaller in the neutral case due to the small temperature 
gradients. Total advection was a cooling process and acted 
in the same direction as diffusion. Due to the smaller 
temperature gradients in the horizontal and the vertical, 
both components of diffusion were smaller in the neutral 
case than in the initially stable one. Slightly faster cool- 
ing resulted overall. 

Differences in the temperature profiles in the 
mid-slope region were small out to 800 seconds or so. 
Slightly larger cooling amounts were found near the ground. 
Aloft, less cooling was predicted in the neutral run as the 


temperature gradients there were quite small. Along the 
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Table 4.6 Contributions to 06/dt from advection and 


diffusion just above mid-slope (J=8,K=7) for Run 6 at 200 
second intervals 
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trough boundary, much reduced cooling was found aloft, 
relative to the standard run, as well, even though the ver- 


tical wind speed was substantially larger. 


4.7.3 More Stable Case Initially - Run 7 
Another model integration, this time under extremely 
Stable conditions, was performed. The lapse rate used 


1 


initially was -0.04 C° m-', creating very large vertical 


potential temperature gradients of about 0.05 C° m™'. In 
this run, however, a counterclockwise cell developed above 
the valley floor as in other runs after only 620 seconds, 
due to the large amounts of cooling there through vertical 
advection. Useful examination of the results can be made out 
as far as this time, at least. The evolution of the minimum 
streamfunction value, shown in Figure 4.14, shows clearly 
that the stability has acted to retard the development of 
the downslope flow (given that the streamline centers were 
located close together). After 600 seconds, for example, the 
maximum horizontal wind speed at the sloping surface was 
O24) e09225 bandrovds'm.st' for’ they runs: with initial lapse 
Pates or 0.01 Co m= (Roneo)., =07004 Come” (Rone 5. and 
-0.04°C° m- 7° (Run 7), °respectively. 

Horizontal advection close to the ground was negative 
and slightly smaller than for the less stable run (Run 1) as 
a result of the somewhat lower wind speeds predicted in Run 
7. Vertical advection, on the other hand, was much greater 


than for the less stable case and warm advection was the 
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overall effect. Changes through diffusion were negative and 
generally somewhat larger than in the standard run due to 
the larger gradients involved. Overall, cooling slightly 
above the surface was smaller in Run 7, mainly as a result 
of the higher rates of warm air advection. 

The variation in: verticalstremperature gradient ewith 
time illustrates the differences in cooling amounts with 
height in different parts of the valley. Table 4.7 gives the 
vertical potential temperature gradient at three different 
valley locations (trough, mid-slope, and ridge) based on the 
temperatures at the ground and 32 m above. Results for each 
of the three initial stability regimes considered in this 
Study are cited. In all cases, the inversion was strongest 
along the ridge boundary where warming by subsidence 
occurred aloft and cooling through diffusion near the ground 
were the dominant factors. Along the trough boundary, cool- 
ing occurred both near the surface and, to a lesser extent, 
further away. This caused some increase in the vertical tem- 
perature gradient with time, although it was not as large as 


that above the ridge. 


4.8 Comparison with Observations - Run 8 and others 

The results presented in this Section are included in 
order to compare and contrast the predicted values and the 
observed. Various initial and boundary conditions were 
utilized, in accordance with the observations. Model 


verification was difficult since the total length of time of 
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Table 4.7 Vertical potential temperature gradients (based on 
surface temperatures and that 32 m above) under different 
lapse conditions for various locations in the valley 
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integration was quite limited. As well, the data base was 
somewhat sparse, in terms of the location and numbers of 
observing stations (typically, three were available - both 
valley rims, plus profiles above some mid-slope region). 
However, some useful comments can be made. 

The observed temperature and wind speed profile 
variations were used to determine the numerical values of 
the many constants involved in the model. For example, the 
evolution of the hourly mean temperatures observed at the 
Dawson Bridge site was used to extract surface cooling 
rates. These screen temperatures, obtained at 1.2 m above 
the ground at the rim and mid-slope regions, are shown in 
Figure 4.15 for three experimental days in 1978. The earlier 
onset of cooling along the slope, compared to the rim, is 
very evident. All seemed to indicate nearly constant cooling 
HaeessoL 2)- 2.oec. “per nour. 

During most obServational evenings, the flux Richardson 
number quickly exceeded its critical limit. As discussed in 
Section 2 of Chapter 2, this brings about total suppression 
of turbulent exchange in the vertical. As a result, (Km), 
was set to zero throughout these final model integrations. 
Initially; while still under nettral conditions, such 
exchange could be significant just as the surface-based 
inversion is forming. In spite of this, a constant value was 
used during the whole period. The vertical exchange of heat 
through turbulence was effectively zero at most times as 


well, However, a non-zero K+ was used here in order to 
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Figure 4.15 Observed evolution of screen temperatures (1.2 m 
above ground) for three experiments in 1978 at the Dawson 
Bridge site at (a.) east rim and (b.) east slope 
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Simulate radiative exchange processes as outlined in Section 
6 of Chapter 2. A value of 0.1 m* s~', estimated by Brunt 
(1934), was used for Kp. 

A preliminary investigation into the effect of 
eliminating vertical diffusion of momentum was made through 
a model integration having all parameters equal to those of 
Run 6, except that (K,,), was zero. Relative differences in 
the minimum streamfunction and maximum surface downslope 
wind speed of about 10 per cent were found after 800 seconds 
of integration time. This was the hoped for result since 
reasonable wind speeds and temperature changes had been 
obtained in Run 6, uSing a vertical exchange coefficient of 
AD Semiees *< O. 

Finally, one integration was performed using the most 
appropriate ‘values of (Kp), = (ky)y = 91/0 m> so, 

(Ky)2 = (Ky). = 0.0, Kp = 0.1 m* s-' “and! ’al constant’ ‘cooling 
rate of 2.5 C° per hour. This run was initiated under 
neutral stability conditions and a valley trough temperature 
of 297.6 °K, and is denoted as Run 8. A counterclockwise 
cell developed high above the valley floor after about 1450 
seconds, as in the other runs. For this reason, the amount 
of cooling predicted after only 1200 seconds (20 minutes) 
was compared with that observed over similar time periods. 
The reduction in the vertical heat transfer coefficient, 
((Ky),2 + Kp), resulted in small amounts of cooling away from 
the surface. At the location of the mid-slope station, a 


temperature profile was obtained on three dates in 1978 at 
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heights of up to 15 m, depending on the particular 
experiment. Based on the hourly mean temperatures, cooling 
amountsS meaSured over the early evening hours from 1900 - 
2100 MDT were computed. Assuming the cooling to be 
distributed uniformly over this period, one can derive the 
approximate amount of cooling in 1200 seconds for comparison 
with the model predictions. Figure 4.16 illustrates these 
observed mean cooling amounts at different heights along 
with that predicted in Run 8 above the mid-slope region. 
Clearly, the agreement is quite good. 

At later times, after 2100 MDT, the amount of cooling 
was observed to decrease at higher levels, perhaps in 
rsponse to effective mixing by the downslope wind. This is 
evident inetne valley modelsresultswwhensthescocling amount 
predicted at some height above the ground Petercrned as a 
function of time. For example, the temperature change 
predicted 8 m above the mid-slope region is 0.015, 0.047, 
0.072, 0.090, 0.096, 0.083, 0.074 C° per 200 seconds based 
on total change over successive 200 second intervals. 
However, integration was not feasible much beyond 1200 sec- 
onds. 

In accordance with the cooling being slowly transferred 
upwards, the downslope flow was confined to a relatively 
thin layer, compared to most other runs. However, due to the 
large horizontal temperature gradients close to the ground, 
considerably larger surface downslope wind speeds were 


predicted (with values up to about 0.9 m s~' after 1200 
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Figure 4.16 Observed cooling amounts versus height in 1200 
seconds (mean values) and that predicted by the valley model 
in Run 8 
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seconds). It is believed that the earlier development of 
Strong winds accounted for the existence of a counterclock- 
wise cell above the trough through large amounts of 
adiabatic cooling there and the positive horizontal tempera- 
ture gradient thus created. 

The profile of the horizontal wind component predicted 
in Run 8 above the mid-slope region is illustrated in Figure 
4.17 (negative values are plotted for flow directed from the 
valley ridge towards the trough). The depth of the downslope 
flow increased slightly as the wind speeds themselves 
developed. Maximum speeds were predicted to occur at the 
ground surface, under the 'free-slip' boundary conditions 
imposed in this model. As well, the flow developed rather 
quickly. Comparison with the mean observed wind shear from 
the valley experiments 1s provided in Figure 4.18. Periods 


during which R, was supercritical were chosen since there 


f 
would be little downward transfer of momentum from a 
prevailing wind aloft during such times. The slope wind was 
measured at a height a 0.8 m uSing a very sensitive Gill 
propeller anemometer. At heights of 1.7, 3.84, and 5.84 m, 
however, Rimco cup anemometers were used to measure the hor- 
izontal wind speed and these had a stall speed of about 0.25 
ms~'. For this reason, only 30 minute intervals with 
average wind speeds above this threshold were used. On the 
two evenings considered, June 27 and July 4, 1978, mean 


maximum speeds of 0.73 and 0.57 m s~', respectively, were 


found at the lowest level of measurement. The predicted 
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Figure 4.17 Horizontal velocity component versus height 
predicted by valley model in Run 8 out to 1200 seconds 
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Figure 4.18 Observed wind shear and that predicted by valley 
model in Run 8 at 400 and 1200 seconds 
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Shear above the slope region, which is very close to that 
predicted at the bottom of the slope where the observations 
were made, agrees quite well with the limited observations 
shown in Figure 4.18. 

A final model integration was performed using (2.6.4), 
the surface cooling rate based on radiational losses. A net 


bediation 19655 “per junitwarea per unitetime; i Ri0t 20 ew man * 


2 - 1 


was used. A vertical heat transfer coefficient of 1.0 m* s 


was implemented in order to include some diffusional and 


some radiational exchange, with Kp = 0.3 m’* 


s”' as 
determined empirically by Anfossi et al. (1976). With these 
and other typical atmospheric values substituted into 
(2.6.4), the surface cooling rate was specified to be 

T(O,t) = T(0,0) - 0.075 t? 
wherewt) is im seconds, and T(zet)) -tsuin. “C..thvs curveris 
plotted in Figure 4.19 out to 1200 seconds, along with the 
cooling produced with the constant cooling rate of 2.5 C*® 
per hour specified in Run 8. Evidently, sizeable cooling is 
specified by (2.6.4) at early integration times although it 
falls off as t°*. In response to this large enforced cooling 
rate initially, rapid development of the circulation 
resulted in downslope wind speeds of about 1.5 m s~' after 
only 600 seconds of integration time. Shortly after this, at 


about 760 seconds, a counterclockwise cell formed above the 


valley trough. 
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T(0,t) = T(0,0) -2.5/3600 t 


~{.0 


T(0,t) = T(0,0) -0.075 t% 


-1.6 


TEMPERATURE CHANGE (C°%) 
-2.0 


0 200 400 soo 48=8si«<iakst=“‘<‘é‘é OC ~C« DOC 
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Figure 4.19 Two boundary conditions used to specify the sur- 
face temperature fields 
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Chapter 5 


SUMMARY AND CONCLUSIONS 


5.1 Summary 

A two-dimensional finite-difference grid model has been 
developed in order to simulate the observed downslope circu- 
lation found to exist within the North Saskatchewan River 
valley. Some simplifications of the equations of motion were 
made, including neglect of Coriolis force, use of K-theory 
to describe turbulent motions, and two-dimensional flow ina 
vertical cross-section of the valley. The assumption of 
incompressible non-divergent flow permitted the use of 
streamfunction fields which are easily related to the valley 
components themselves. The thermodynamic equation, coupled 
with the vorticity equation, was solved numerically to 
predict the velocity and temperature at each grid point. 

Downslope winds are produced in direct response to hor- 
izontal pressure gradients created through specified cooling 
rates at the ground. Radiative heat transfer was 
approximated in the thermodynamic equation as a diffusion 


process, with a radiative diffusivity, K This followed 


R° 
considerations involving long-wave absorption and 
re-emission by water vapour and was introduced by Brunt in 
the 1930's. Two means of stipulating the surface boundary 
conditions were used. First, constant surface cooling rates 


were simply applied, as observed within the valley under con- 


Sideration. Secondly, the radiative heat transfer equation 
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was solved analytically to predict an equivalent rate of 
temperature change at the ground. As this is the driving force 
for the development of the valley circulation, the exact nature 
of the temporal variation of the ground temperature is most 
Significant. 

Much of the model development in this study dealt 
with the numerical aspects of the finite-difference 
formulations of the terms in the governing equations. 

Due to problems inherent in the Dufort-Frankel scheme, 

when applied to diffusion terms in equations possess- 

ing several types of terms, its usSe was finally abond- 

oned in favour of a forward-time, centered-space 
formulation. Advective terms were handled by a forward-time, 
upstream-space configuration. Solution of the elliptic 
partial differential equation relating streamfunction to 
vorticity was obtained using Liebmann sequential relaxation. 
This procedure consumed most of the CPU requirements. 

Each integration was initiated with the atmosphere at 
rest .wAS avresudit!( turbulentedifiusione and’ radiation 
processes dominated in the early stages of each run. As the 
Circulation developed, advection became increasingly 
important, depending upon the particular application. For 


* s~'predicted 


example, one run with (K,), = (Ky) 2 =n. 00m 
rates of change of temperature by advection and diffusion/rad- 
lation ustmabovel thetislope ofr sseixia 0 geandt=-8 xf 10h’ ceys;, 

respectively, whereas another with (K,,), = (Ky), = 0.01 m* s7' 
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integration times (1200 seconds). This illustrated how ad- 
vective and diffusive processes both strive to destroy any 
temperature anomalies present. When diffusive heat transfer 
is necessarily limited by the presence of marked stability, 
advective processes (whose magnitudes depend partially on 
gradients present) are proportionally increased. 

The prediction of, or at least hint at the eventual 
development of, a counterclockwise flow high above the flat 
valley floor was evident in each integration. The timing of 
1tS appearance varied from run to run, however. This cell 
resulted from differences in vertical velocity (and thereby 
amounts of cooling through adiabatic expansion) at the same 
height across the valley, with a cooler area being produced 
close to the trough boundary. Through the solenoid term in 
BHeeVOrtrcity requation (which Ys “proportional to *the hori 
zontal gradient in potential temperature), increasing vorti- 
city was produced, and eventually positive vorticity 
values were attained in this region. A counterclockwise 
cell was the result. Under conditions of extreme 
Stability, this phenomenom occurred quite early (after 
about 600 seconds), i.e. even very low vertical wind 
speeds combined with large vertical temperature gradients 
could produce a differential in cooling amounts at the same 
height at small integration times. 

The sensitivity to the particular geometry of the val- 
ley cross-section being modelled was examined. In most runs, 


the upper boundary was located at a height of 108 m, close 
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to twice the valley depth. However, one integration with a 
lid at 144 m suggested that the evolution of the thermody- 
namic and flow fields was not very much affected, at least 
within the depth of the valley. As well, most integrations 
were performed within a cross-section which featured a flat 
valley bottom. When excluded, the flow was forced to rise 
sharply at the base of the slope and relatively higher ver- 
tical wind speeds resulted. In this case, the counterclock- 
wise cell developed somewhat earlier. The conditions near 

the sloping surface and the rim were largely unaffected. One 
integration was executed with the sloping surface at an 

angle of 0.197 radians (compared with a usual value of 0.100 
radians). In this case, downslope winds were up to 70 per 
cent higher, depending upon time of comparison. However, this 
difference was not maintained and, after 800 seconds, maximum 


-' were predicted along the 


surface slope winds of 0.37 ms 

steeper slope whereas in the standard run, they were 0.34 ms-“'. 
The magnitudes of the downslope wind speeds under a 

variety of initial stability regimes were investigated. As 

in other studies (both observational and modelling types), 

the development of the valley circulation was severely 

hindered in very stable situations. Detailed comparison with 

experimental data obtained within the North Saskatchewan 

River valley was difficult but generally good agreement was 

found in the low-level temperature profiles above the slope 


and in the magnitudes of the slope winds themselves. The 


exact time scale for the evolution of the valley wind and 
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temperature regimes is hard to ascertain. The observed wind 
reversal (decrease in speed with height) was predicted ; the 
model (under conditions of 'free-slip' at the lower boun- 
dary) produced maximum wind speeds coincident with the 
ground surface. Increased stability was brought about at the 
rim poSition, relative to the valley bottom, through 


subsidence. 


5.2 Conclusions 

The numerical model described some observed propert- 
ies of the micrometeorology of the North Saskatchewan 
River valley. Although total integration times were 
only 10 - 20 minutes, useful conclusions regarding the 
development of the downslope flow can be made. Integ- 
Pipi ens were carried out by Thyer (1966) using a similar 
finite-difference representation of the equations over 
a Simple 'V-shaped' region out to only 120 seconds. 
The advent of computational instability difficulties 
may have precluded longer integrations although this was 
NOt Statea Explicitly by therauthor.. The= inclusion of a 
flat valley bottom into the present model provides a 
more realistic simulation for applications to the North 
Saskatchewan River valley. 

Investigations into the various possible numerical 
formulations of the terms in the governing equations led to 
the choice of forward-time, upstream-space for the advection 


terms and centered-time, centered-space for the solenoid 
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term in the vorticity equation. The appropriate selection of 
a scheme to approximate the diffusion terms was more 
Gifficult ; however, after applying several schemes to a 
Simple one-dimensional diffusion equation, for which the 
analytical solution was derived, a two-level scheme in time 
proved desirable in this application. This resulted from the 
presence of other terms in the equations, which caused a 
disturbance in the computation of diffusive changes through 
altering the values at intermediate times. 

Through variations in initial stabilities, surface 
cooling rates, and vertical transfer coefficients, the 
relative .ccontmibutions -of advect hon, “aiifusionprand 
radiation could be estimated under several regimes. Results 
from one integration, which was meant to simulate observed 
valley inversion conditions, clearly demonstrated the 
importance of radiative cooling processes close to the 
ground in the early stages of the slope wind development. At 
a location 4 m above the mid-slope region, for instance, 
radiational cooling was at least an order of magnitude 
larger than advectional cooling, at 200 seconds of integra- 
tion time. Gradually, as the slope winds increased, changes 
in temperature due to advection increased. Horizontal 
advection acted to transport air away from the cooling 
surface, whereas vertical advection brought down warmer air 
from aloft. At all times, the net result was negative tem- 
perature advection although the differences in the 


magnitudes of these two terms were small. Eventually, after 
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1600 seconds of integration time, advective changes were 
about twice as large (-4 x 10°“ C° s°-') as radiative or 
diffusive ones. 

For comparison, the analytical solution to the one-di- 
mensional vertical diffusion equation, using a diffusivity 
of 0.1 m* s~' equal to that used in the above valley model 
integration, predicted a very similar low-level temperature 
profile. Therefore, when advection processes were not 
included, larger gradients resulted in increased radiative 
cooling and it seems evident that radiative heat transfer 
alone was effective in this application. 

The model predicted quasi-steady-state conditions 
within about 20 minutes of integration time with 
development of downslope wind speeds of = 0.8 ms~'. Both 
the predicted low-level temperature variations with height 
and the maximum slope winds agreed extremely well with the 
limited observations made within the North Saskathewan River 


valley in Edmonton. 


5.3 Suggestions for Future Work 

Due to the large amounts of computer run-time and 
storage costs, the number of model integrations was 
obviously restricted. Experimentation with various smaller 
and larger grid sizes was not feasible, and this may cast 
some uncertainty on the validity of the results. As the 
relaxation procedure appeared to consume the majority of the 


total cost for each integration, investigations into other 
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numerical or analytical methods of solving elliptic partial 
differential equations would prove extremely valuable. Even 
the Liebmann sequential relaxation technique used possesses 
various procedures by which the iteration process can be 
speeded up. For example, Haltiner (1971) indicated that 
convergence was more rapid over a N. Hemisphere grid when 
the sequence of points followed a pattern which began at the 
outer boundary and spiralled inwards. Similar strategies may 
have proven useful in this application. Also, the optimal 
over-relaxation coefficient may have varied from run to run. 

The use of Fickian diffusion (i.e. constant in space 
and time) may be a major drawback in this study. The 
apenOpriate values Of Ky sand Keane intimately related to 
the variations in the wind and temperature fields, and these 
continually evolve as the integration proceeds. Coefficients 
which depend on the local stability would be important in 
this application. In addition, transfer properties close to 
the lower surface are likely to be much different from those 
well above the ground. The incorporation of radiative 
transfer was initially considered to be beyond the scope of 
the study, but a simple approximation by a diffusion process 
permitted its ultimate inclusion. 

Provision of an area of integration consisting of both 
slopes which define a given valley, using lower boundary 
conditions which simulate earlier shading on one Secs than 
on the opposite one, would provide insight into the 


asymmetrical nature of the resultant flow. The 
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two-dimensional treatment of the problem may also be 
Significant if the circulation predicted in this study were 
to be used as input into a pollutant transport model such as 
Rudolph (1980). In this case, a downvalley wind would create 
a helical trajectory with recirculation of pollutants 


possible if the source extended upvalley some distance. 
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APPENDIX A ;: STABILITY ANALYSIS FOR ADVECTION EQUATION 


The forward-time, upstream-Sspace finite-difference 
scheme is examined for stability when applied to a one-di- 


mensional linear advection equation such as 


3S 3s 
ae Yinean 42° ' | CA 


The von Neumann or Fourier Series method is most often used, 
according to Mesinger and Arakawa (1976). Letting 
S(y=JAy,t=ndt) = S(J,n), (A.1) can be written using 
finite-differences as 

s(J;nti) -<S(J,n) _ S(J,n) - S(J-1,n) 

At - ~v | Ay eee : 

or SCisutt)s = (1-1) S(i,n) +) SQ—15n)) 4) v0> 70 ’ (A.2) 
where uw=v At/Ay. 


Substituting a solution of the form 


S(J,n) = S(0,n) exp(igJAy) , (Arse) 

one obtains S(0,n+l) = (1-y) S(0,n) + wp S(0,n) exp(-ifAy) « (A.4) 
An amplification factor, A, 1s defined such that 

S(O,ntl) = |A| $(0,n) ‘ (A.5) 


The von Neumann criterion for stability can be stated as 


cee (A.6) 


Substitution of (A.5) into (A.4), and simplification, yields 


A = d—t ten,exp(-LBAy) 


143 


a ee " Vin 
elem, 
“WOrYA0OS 


a 


wor" 
enne19t2lb etial? _v. — ovtequ « 
-~z&-spno s 07 56 —— nee erations i a9 
7 >. 


a 


as Se ups ag i be 


fi ae 
_ 
2 
e er ia 
; bom & sonljem @eties Wei ww 
‘ Ewasé iG 


, ely ’ 


ea euins enti ibs 


Tea en eee wee a Citaalye 


= 
to nolawiog 90 oie 


+ ee al 


(valgslexs (e,b)2 « in We a 

Mit-)luxs (o,0)2 ¢ + (a,@)8 Ewel) © (1e0,0)8 teide 
bentiah ak ,4 ,2090e2 motze s3itkiqna i 

; (a,0)2 [A] = Che Ope 7 


. . 
: badese ¢ on Yiltdeda tol eeliszite en 


‘ te i 


(yAG be) inne F Yi 


144 


Jal? = (1 - up + up cos(BAy))2 + (-n sin(Bdy)) 


[al@ = 1 - 2n(1 - ») (1 =- cos(BAy)) ; GAN) 


Thus, von Neumann's criterion for stability becomes, from 


(A.6), 


1 = 2n(1 - up) (1 - cos(BAy)) < 1 : (A.8&) 


Since u > 0 and i-cos(BAy) > 0, (A.8) becomes simply 
yFv—<l : (A.9) 


This is the well-known Courant-Friedrichs-Lewy (C.F.L.) 
Eriterion. 

The von Neumann method of examining stability based on 
Fourier Series, is valid only if the coefficient of the 
difference equation (v in (A.1)) is constant. The method can 
be applied locally if the difference equation has a variable 
coefficient. Mitchell and Griffiths (1980) cite that there 
is much numerical evidence to support the contention that a 
method will be stable if the von Neumann condition is 
Satisfied at every point of the field, even if derived as 
though the coefficients were constant. 

To obtain further information on the behavior of the 
numerical solution, (A.8) can be examined in more detail. A 
sketch of |A|? versus uw for various wave numbers is given in 


Figure A.1. Within the stable region (|A|* < 1), this scheme 
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is damping for all values of uw < 1 and the amount of damping 
depends on the wavelength. In this derivation, however, the 
true solution has a constant amplitude - this indicates some 
type of error due to the finite-difference representation. 
These problems should be kept in mind when the results in 


Chapter 4 are considered. 
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Figure A.1 Amplification factor versus uw using forward-time, 
upstream-space scheme applied to advection equation 


ills ,3evawort ,nol 
amos moore te ‘yi ob 


-foitsiaseeige? po i 


nl edices: ot? “nedy om 


yas\e@ 8 


.amjg~biswro! entaw « euerev togae2 nolzasiiite 
noljazpe ootscevhe os Gebiggs saaioa 4 


APPENDIX B - STABILITY ANALYSIS FOR DIFFUSION EQUATION 


The Dufort-Frankel finite-difference scheme is first 
examined for stability when applied to a one-dimensional 


linear diffusion equation such as 


9S _K 4S 
rye 5 oy> e CBE) 


The von Neumann method of stability analysis is applied, 
uSing the same notation as in Appendix A, as given in an 
exercise in Mitchell and Griffiths (1980). (B.1) can be 
written in finite-difference form as 

S(J,nt+1) (1 + 2r) = 2r {S(J+l,n) + S(J-1,n)} + (1 = 2r) S(J,n-1) (B.2) 
where r = K At/Ay*. Rewriting this as a two-level system, in 


vector form, one obtains 


S(J,n+1) (2r/(1+2r)) D (1-2r) /(1+2r) S(J,n) ee 
- Be 
T(J,nt1) ub 0 ECs, 1) 
where D = centered difference operator such that 
D S(QWJ,n) = S(J+1,n) - S(WJ-1,n) 
Substituting a solution of the form 
U(J,n) = U(O,n) exp(iBJAy) ’ (B.4) 
where 
S 
U = 
As 
and rewriting the result gives 
(4r/(1+2r)) cos (BAy) (1-2r)/(1+2r) 
U(O ntl) = U(O,n) (B.5) 
a 0 


amplification matrix 
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The elgenvalues, A;,, of the amplification matrix above 
Satisfy |A - AI| = 0 and are found to be such that 


1.2 - 4x cos(BAy) A, - (1-2r) 
de es eet Aa 


1+2r 1+2r 


Solving this quadratic in A;, one obtains 


2r cos(BAy) + (-4r2 sin (BAy) +1)? 


17 (B.6) 


The von Neumann necessary condition for stability of a 
two-level system is 
max |A,| <1 (B.7) 
From (B.6), the maximum value of A, or Az iS seen, by 
mispection ms tossatisty (8.7) forall valuesr ofrm="ke ary Ay. 
Thus, the Dufort-Frankel scheme is unconditionally stable. 
The forward-time, centered-space scheme will be 
examined in terms of stability also using the von Neumann 
method. This scheme, applied to (B.1), can be expressed as 
S(J,n+1) = (1-2r) S(J,n) + r{S(J-1,n) + S(J+1,n)} GBs) 
where r = K At/Ay”? as before. Substitution of a solution of 
the form 
S(J,n) = S(0,n) exp(ifJAy) OX ey, 
results in, after simplification, 
S(0,n+l) = S(0,n) {(1-2r) + r exp(-iBAy) + exp(iBdAy) } (A 24) 
Using the definition of an amplification factor given by 
(A.5), yields A} = 1 - 2r + rfexp(-iBAy) + exp(ifBdAy) } 


According to (A.6), numerical stability results in 


[Al= [1 - 4r sin*(BAy/2)| < 1 
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By inspection, this occurs with 0 s r s 1/(2 sin? (BdAy)) 
which leads to the stability condition 


K At 


ai 
— (Be 9) 
Ay? Z 


((yst) "ate SINE 
fee 


(2.8) 


APPENDIX C - ANALYTICAL SOLUTION TO ONE-DIMENSIONAL 


DIFFUSION EQUATION 


In this appendix, an analytical solution to the one-di- 


mensional heat diffusion equation defined by 


30 926 | 
—=K—,0<2z2<H,t>Q0 f G35) 
ot 922 


is derived, subject to the following boundary conditions. 


e(H,t)=T,t>oO , (355 12) 


6(0,t) = 6(0,0) - at, t > 0 : GER ace, 


and initial conditions 


Oz, O)m= (Cr. —10:(0;0))oz7Hi+"6(0,0).).0 <22eeen oe (35659 
Following the methods for solving partial differential 
Seiatenco risierieh non-homogeneous boundary conditions suggested 
in Boyce and di Prima (1965), the following transformation 
is made : A(z,t))= 0(z, C)o—sb Czy) Comets) 
Therditferential. equation, (3.5.4) aun terms Of A and: Bis 
now 
dA dB 397A 9°B 
—+—-K—-K—=0,0ez2<«H, t>0 : (Gz.2)) 
dt at 9z* 92° 

and: tthe *boundaryrconditi onsiee(Si5 52 )irandie (3.2573) - famennow 


A(HG.t). +) BCH, t).. =*T,. tr>s0 ‘ CEA) 


A(O,t) + B(O,t) = 6(0,0) - at, t > 0 : erg) 
The initial conditions, (3.5.5), become under this 


transformation 
A(z,0) + B(z,0) = (T — 6(0,0))z/H + 6(0,0), 0O< z<H ee Ome 
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It is now required to find some B(z,t) such that B(H,t) 
= T and B(0,t) = 6@(0,0) - at in order that the boundary con- 
ditions, (C.3) and (C.4), become homogeneous. By inspection, 


this can be accomplished if we set 


B(z,t) = (6(0,0) - at) (1 - z/H) + zT/H ° (CAG?) 


Withethisechoice,of Biz, fo. 


22 = - a(1 - 2/H) 5 
9°B 
K —— = 0 2 
922 


and so (C.2) becomes 
oA 072A 
— - K — =a(l1 - z/H) = q(z), O< z< H, t>0 : CGT.) 
ot az 

and as well, the boundary and initial conditions ((C.3), 


(c.4), and (C.5), respectively) become 


ACHet)E=105, ts sO ae (Cx 8}) 
A(Q,t).= 0, t.> 0 ’ (C.9) 
A(z,0) = 0, O< z<H : (Cou) 


Note that the above boundary conditions are now homogeneous 
as desired. However, the differential equation, (C.7), has 
become non-homogeneous. A second transformation is therefore 
required : C(z,t) = A(z,t) - D(z) i Coalalt) 
The differential equation, (C.2), can now be written as 


9c 92C 92D 
ee (ie) j (Ce 12) 
ot gz? 92° 


By setting 
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92D 
-~K—=q(z), O< z<¢H, t > 0 : (C213) 
gz- 
equation (C.12) becomes homogeneous as required, with the 
corresponding boundary conditions 


D(H) = 0 9 (Ema) 


D(O) = 0 ‘ Gra oo) 


From (Cc. 7), g(z)"="e W1e— 2/H), So (Cl13) becomes 
22D Oa 
— es-(z/H-1),0<z2<H : CC al col 
gz K 
The general solution to this, subject to (C.14) and (C.15), 
is 
a 
D(z) = — {z9/6H - z*/2 + Bz/3}.. (Ga) 
K 
found by integrating (C.16) and evaluating the constants of 
integration using the boundary conditions. Substituting this 
back into (C.12) results in a homogeneous equation in C(z,t) 
alone ;: IC 92C 
—-K—=0,0<2z2<¢H,t>0 ' Ger t.G.) 
with boundary conditions of 


C(Het) =a0e te>:.0 ; CC. 19%) 


ei Cipaeyy = 18) Asa dv , Geen 


Ang initial conditions of 


C(z,0) = A(z,0) - D(z) = - 3 {z3/6H-z2/2+zH/3} = p(z) CGa2)) 
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The general solution using separation of variables 
technique is found in any partial differential equations 
textbook, to be 


C(z,t) = } b| exp(-n?n?Kt/H?) sin(nnz/H) , O< z<H, t>0 (C.22) 


n=] 


where b 


2/H frp (z) sin(nmz/H) dz 


a 
2/H [- — {23/68 - 22/2 - 2H/3} dz 
K 


= - 2H2q/ (nm) 3K 
Backtransforming twice to return to the original variable, 
ONG Ze Ge SOG oc) eC lait) e im Ci ciste) Pern C Ce) enra Cw Gaus) ie CV.e,S 


the general solution to (3.5.1) to be 


fos) 


@(z,t) = ) { (-2H2a/ (n7) 3K) exp(-n¢1*Kt/H2) sin(nmz/H} 
n=1 


a 
ti-aw(z9/6Hezs/2 -azH/3) +-(6(080))— at) (1 =--27H). (3.52.6) 
K 


+ 2T/H,; Ons Ace, .t > 0 
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APPENDIX D - SOURCE CODE FOR COMPUTER PROGRAM 


The following source code was developed by the author 
for execution on the University of Alberta's AMDAHL 580/5860 


computer, which operates under the Michigan Terminal System. 


SS SL TE a SE Sa ESS TE DOD ENE IR GE FP EE TE TS TE SITE EE ET I TEE TEE 


FMPETCTT PREAU*S CAS H..0 == 2) 

DIMENSION THO(16.28), ZO(16,28), VO(16,28), WO(16,28), THi( 16,28), 
210 16,28), Vi01G,28), wi€i6,28), PHI( 16,28), A1(16,28). 
A2(16,28), T(28), ZKM(16), TKM(16), A3( 16,28), 
A4(16,28), A5(16,28), A6(16,28), FOR1(10), FOR2(10), 
ZKM2(16), TKM2(16) 

INTEGER TEND 

REAL*8 MAGNV, KHZ, KVZ, KHT, KVT 

COMMON TTOT, FORi, FOR2, KL, JEND, JCOR, KEND, KMID, JJB 

COMMON /A/ KHZ, KVZ, KHT, KVT 


AWN = 


Cc 

Cc THE FOLLOWING ARRAYS ARE DEFINED OVER THE VALLEY CROSS-SECTION 

Cc THO - POTENTIAL TEMPERATURE AT PREVIOUS TIME STEP 

Cc TH1 - POTENTIAL TMPERATURE AT PRESENT TIME STEP 

Cc ZO - VORTICITY AT PREVIOUS TIME STEP 

Cc Z1 - VORTICITY AT PRESENT TIME STEP 

Cc PHI - STREAMFUNCTION AT PRESENT TIME STEP 

Cc VO ~- CROSS-VALLEY WIND SPEED AT PREVIOUS TIME STEP 

Cc V1 - CROSS-VALLEY WIND SPEED AT PRESENT TIME STEP 

Cc WO - VERTICAL WIND SPEED AT PREVIOUS TIME STEP 

C Wi - VERTICAL WIND SPEED AT PRESENT TIME STEP 

Cc Ai - ADVECTIVE CHANGES TO POTENTIAL TEMPERATURE 

Cc A2 - DIFFUSIVE CHANGES TO POTENTIAL TEMPERATURE 

Cc A4 - ADVECTIVE CHANGES TO VORTICITY 

C AS - DIFFUSIVE CHANGES TO VORTICITY 

C A6 - SOLENOID TERM 

Cc LOGICAL UNIT 5 : INPUT 

C LOGICAL UNIT 6 : POTENTIAL TEMPERATURE, A1,A2 

Cc LOGICAL UNIT 7 : STORES LAST TWO TIME STEPS FOR CONTINUATION OF 

C RUN AT A LATER TIME 

C LOGICAL UNIT 8 : VORTICITY, A4,A5,A6 

Cc LOGICAL UNIT 9 : STREAMFUNCTION 

Cc LOGICAL UNIT 10 : CROSS-VALLEY WIND SPEED 

Cc LOGICAL UNIT 11 : VERTICAL WIND SPEED 

C LOGICAL UNIT 12 : DATA FOR LAST TWO TIME STEPS WHEN RUN BEING 

C CONTINUED AT A LATER TIME 

C 
CALL "FPREAD(US? “1G4R*872T:", TEND, EPS. RAEPHAY TC1). OELTAT? KEND; 
1 MAXITR) 

Cc 

Cc TEND - FINAL INTEGRATION TIME DESIRED 

C EPS - CRITERION USED TO STOP RELAXATION PROCESSS 

C “RALPHA - OVER-RELAXATION COEFFICIENT 

Cc T(1) - INITIAL SURFACE TEMPERATURE 

Cc DELWAIe — EENGTHSORS TIMES STEP 

Cc KEND - HEIGHT OF UPPER BOUNDARY (GRID UNITS) 

Cc MAXITR - MAXIMUM NUMBER OF RELAXATION PASSES ALLOWED 

Cc 


CALL FREAD(S ~ (ar as. 0217, D022, DY) 


DZ1 - VERTICAL GRID SIZE IN LOWER PART OF VALLEY CROSS-SECTION 
0227 —" VERTICAL. GRIDVSIZE] IN] UPPER PART OF SVAELEY CROSS-SECTION 
DY - HORIZONTAL GRID SIZE EVERYWHERE 
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RATIO = DZ2 / D0Zi 

IRATIO = RATIO 

KVLY = (52. + DZ1)7/ 021 
KVLYMi = KVLY - 1 

KMEDo =) C72. *202Z1) 7/7024 
JEND = (600. + DY) / DY 
KENOMI = KEND - IRATIO 
JCOR = (80 + DY) / DY 


KVLY - HEIGHT OF RIDGE (GRID UNITS) 

KMID - HEIGHT AT WHICH VERTICAL GRID SIZE IS REDUCED (GRID UNITS) 
JEND - LENGTH OF HALF VALLEY CROSS-SECTION (GRID UNITS) 

JCOR - LENGTH OF FLAT VALLEY FLOOR (GRID UNITS) 


DATA G /9.806D0/, CP /1012.0D0/, RHO /1.2D0/, 


ALPHA = DATAN(O. 1D0) 
CALL FREAD(5, ’2R:’, GAMMA, 


COOLR) 


THE FOLLOWING ARE PHYSICAL CONSTANTS 
G - ACCELERATION DUE TO GRAVITY (M/S/S) 


CP - SPECIFIC HEAT OF AIR AT CONSTANT PRESSURE (JU/KG/DEG C) 


RHO - DENSITY OF AIR (KG/M/M/M) 


R - GAS CONSTANT FOR DRY AIR (JU/KG/DEG C) 

ALPHA - ANGLE OF SLOPING GROUND SURFACE (RADIANS) 
GAMMA - LAPSE RATE OF TEMPERATURE (C DEG/M) 

COOLR - RATE OF COOLING AT SURFACE (C DEG/S) 


CALLA FREEADS. e021 3° ,.°K 
READ (5.790) (FORT (1). 
READ (5,790) (FOR2(1I), 
CALE. FREAD GS. 1). 2R*8:; 
ITPR = TPR 


ee 
= uh 


KL, JJB - DETERMINE PORTION OF GRID TO BE PRINTED OUT BY S/R 


OUTPUT 


10 


20 


30 CALL SUBK(JEND, KEND, JCOR, 


KOPT - IF O : DIFFUSIVITIES ARE CONSTANT EVERYWHERE 
IF 1 : DIFFUSIVITIES CAN BE DEFINED 
USING S/R SUBK(JEND,KEND,JCOR,KVLY) 
TPR - TIME INTERVAL AT WHICH MATRICES ARE PRINTED OUT 


JUB) 
10) 
“AOD 
KOPT, TPR, 


COPT ) 


COP, = FO) seINTEGRATIONSBEGINS Alen) = 


CALL FREAD(S, ‘’4R*8:’, KHZ, 


KVZ, KHT, 


KHZ - HORIZONTAL MOMENTUM DIFFUSIVITY 


KVZ - VERTICAL MOMENTUM DIFFUSIVITY 


KHT - HORIZONTAL THERMAL DIFFUSIVITY 


KVT - VERTICAL THERMAL DIFFUSIVITY 


DO 10 K = 1, KEND 
DO .j0).U..2 JEND 

VO(J,K) 0.DO 
WO(JU,K) 0.DO 
ZO(UJU,K) o.DO 
PHI(JU,K) 0.DO 
Vi(JU,K) 0.DO 
Wi(JU,K) 0.DO 
ZACS KD) 0.DO 

CONTINUE 

DO.20. J = aly UEND 

ZKM(JU) = 0.DO 

IF (KOPT .EQ. 1) GO TO 30 


Yow wh 


ITC =O 

IP1 =.0 
iON .=4O. DO 
GO TO 40 


KVLY) 


R /287.0D0/ 


O SECONDS 
IF 1 : INTEGRATION CONTINUES ACCORDING TO LOGICAL UNIT 


KVT) 
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70 


80 
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100 
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130 


140 


TEMPERATURE ABOVE TROUGH (INITIALLY) 


IF (COPT .E€Q. 1) GO TO 130 

DO 70 K 2, KEND 

IF (K .GT. KMID 
GO TO 70 


.AND. (K 


IF (K .GT. KMID) GO TO 50 


INC = 1 

GO TO 60 

INC IRATIO 

T(K) LtKke- 
CONT INUE 


INC) 


POTENTIAL TEMPERATURE ABOVE TROUGH - 


C1 = G / CP - GAMMA 
THO( 1,1) TCT) * 
222 O.D0 
P 1000. 
THO( 1,2) 
Z D21 
P 1000. -* -(TC2)/THO(1, 2) 
DO 100 K = 3, KEND 
IF (K .GT. KMID 
GO TO 100 
IF (K .GT. 
TNC e204 
DZ DZ1 
GO TO 90 
INC IRATIO 
BZ 27022 
THO( 1 ,.K) 
Z, (K*DZ1) 
P 1000. 
CONTINUE 


EME CT) /THOCTS 1) 


.AND. (K 


= 18)74 4) 


) 


) 


KMID) GO TO 80 


1 


(1000. /930. ) 


xk & 


DEXP(DLOG(THO(1,1)) 


* 


4 


DEXP(DLOG(THO(1,K - 


4° CT OK), THOC1:,.K) ) 


)/RATIO - (K - 


- GAMMA * DZ1i * INC 


INITIALLY 


#¥ a (R/ CR) 


(CP/R) 


1)/IRATIO .NE. 


sis, 


0.0) 


e C1024 COTO Wit 1027) 427)) 


(CP/R) 


)/RATIO - (K - 


UNG =2)) 242 C120 Z* 20 TKa— 


Pe CCP AR) 


POTENTIAL TEMPERATURE WITHIN VALLEY CROSS-SECTION 


DO 120 K ip 
JFINAL K - 
IF (JFINAL .GE. 
DO 110 J = 2,. JFINAL 

IF (K .GT. KMID .AND. 
GO TO 120 
THO(UJU,K) THO COs — 
CONTINUE 

CONTINUE 

IUNIT = 6 

WRITE (IUNIT,690) TTOT 

CALL OUTPUT(IUNIT, THO) 

IUNIT 8 

WRITE (IUNIT,730) TTOT 

CALL OUTPUT(IUNIT, ZO) 

GO TO 160 

IUNIT = 12 

READ (IUNIT,850) ITC, 

TF RCL. LO eal) (Goer ON 150 

READ (IUNIT,140) IDUM 

FORMAT (11) 

CALL. UNSAVE(IUNIT, THO) 

READ (IUNIT,140) IDUM 

CALL UNSAVE(IUNIT, ZO) 

READ (IUNIT,140) IDUM 

CALL UNSAVE(IUNIT, PHI) 

READ (IUNIT,140) IDUM 

CALL UNSAVE(IUNIT, VO) 

READ (IUNIT,140) IDUM 

CALL UNSAVE(IUNIT, WO) 

GO TO 160 


1 + JCOR 


1,K 


(K 


) 


ILO 


JEND) JFINAL 


JEND 
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150 READ 


160 


170 


180 


(IUNIT,140) IDUM 
UNSAVE(IUNIT, THO) 
(IUNIT,140) IOUM 
UNSAVE(IUNIT, TH1) 
(IUNIT,140) IDUM 
UNSAVE(IUNIT, ZO) 
(IUNIT,140) IDUM 
UNSAVE(IUNIT, 21) 
(IUNIT, 140) IDUM 
UNSAVE(IUNIT, PHI) 
(IUNIT,140) IDUM 
UNSAVE(IUNIT, VO) 
(IUNIT,140) IDUM 
UNSAVE(IUNIT, V1) 
(IUNIT,140) IDUM 
UNSAVE(IUNIT, WO) 
(IUNIT,140) IDUM 
UNSAVE(IUNIT, W1) 


CALL 
READ 
CALL 
READ 
CALL 
READ 
CALL 
READ 
CALL 
READ 
CALL 
READ 
CALL 
READ 
CALL 
READ 
CALL 


CHANGES IN POTENTIAL TEMPERATURE AND 


TTOT 
ITTOT 
DO 320 K = 1, KEND 
IF (K .GT. KMID 
GO TO 320 
JFINAL K + JCOR - 1 
IF (UFINAL .GT. JEND) JFINAL 
IF (K .GE. KMID) GO TO 170 
DZ DZ1 
INC = 1 
GO TO 180 
DZ = DZ2 
INC IRATIO 
DO 310Ud, = 1, 
IF (K 
IF (K 
IF (K 
IF (K 
IF (J 


TTOT + DELTAT 
TTOT 


.AND. (K - 1)/RATIO 


JEND 


JFINAL 
KMID) ZKM(J) 
: KMID) TKM(U) 
s{a@hs ae olekee Ie slstels ow) 
: KEND) GO TO 220 
SEQ] ap WOR. W)). EO. JEND) ©GG 


Z2KM2(J) 
TKM2(J) 
- JCOR 


INTERIOR 


LES (LTC S60) 4) GO. 10.490 
DTHDY = (THO(JU + 1,K) - THO(UJ 
CALL ADVEC(ZO(JU + 1,K), ZO(UJ,K)., 
HADV ) 
CALL ADVEC(ZO(U,K + INC), 
VADV) 
A4(JU,K) = 
DIFY = 
KHZ 
DIFZ = (ZO(JU.K + INC) 
DZ * KVZ 
A5(JU,K) 
A6(U,K) 
Zi(JU.K) 
ZKM(JU) = ZO(U,K) 
IF (K .EQ. KMID - 
CALL ADVEC(THO(uU + 
HADV ) 
CALL ADVEC(THO(U,K + INC), 
WO(JU,K), VADV) 
Ai(Ju,K) -(HADV + VADV) 
DIFY = (THO(JU - 1,K) - 2.DO*THO(UJ, 
DY * KHT 
DIFZ = (THO(U,K - 
f- OZ 8 KV 
A2(JU,K) = 
THi(JU,K) = 
TKM(JU) = THO(J,kK) 
IF (K .EQ. KMID - 
GO TO 310 


1, 
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-(HADV + VADV) 
C200 se = ai Ky) 


IBN Se MBNA 
G / THO(U,K) * DTHDY 
ZO(JU,K) + DELTAT * 


“wom ow 


IRATIO) ZKM2(uU) 
17K). THOU Ks 
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VORTICITY 


- (K - 1)/IRATIO .NE. 0.0) 


th) GOW 10, 200 
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INTERIOR OF REGION AWAY FROM BOUNDARIES - SUBSEQUENT TIMES 


AT 


AT 


DTHDY = (TH1i(U + 1,K) - THO(U - 1,K)) / DY / 2.00 

CALL DIFFOR(ZO(JU + 1,K), ZO(JU,K), ZTEMP, DY, KHZ, DIFY) 

CALL DIFFOR(ZO(JU,K + INC), ZO(U,K), ZKM(J), DZ, KVZ, DIFZ) 
A5(JU,K) = DIFY + DIFZ 

A6(J,K) = G / THi(JU,K) * DTHDY 

CALL ADVEC(ZO(JU + 1,K), ZO(JU,K), ZTEMP, DY, VO(U,K), HADV) 
CALL ADVEC(ZO(JU,K + INC), ZO(JU,K), ZKM(JU), DZ, WO(JU,K), VADV) 


A4(J,K) = -(HADV + VADV) 

ZTEMP = ZO(J,K) 

ZO(JU,K) = Z1(JU,K) 

Z1(U,K) = ZTEMP + 2.DO * DELTAT * (A4(JU,K) + A5(JU,K) + AG(U,K) 


) 

ZKM(JU) = ZTEMP 

IF (K .EQ. KMID - IRATIO) ZKM2(JU) = ZTEMP 

CALL ADVEC(THO(JU + 1,K), THO(JU,K), THTEMP, DY, VO(U,K), HADV) 

CALL ADVEC(THO(JU,K + INC), THO(U,K), TKM(JU), DZ, WO(U.K), 
VADV) 

At(J,K) = -(HADV + VADV) 

CALL DIFFOR(THO(U + 1,K), THO(JU,K), THTEMP, DY, KHT, DIFY) 

CALL DIFFOR(THO(JU,K + INC), THO(JU,K), TKM(JU), DZ, KVT, DIFZ) 

A2(JU.K) = DIFY + DIFZ 

THTEMP = THO(J,K) 

THO(JU,K) = TH1i(JU,K) 

TH1(JU,K) = THTEMP + 2.D0 * DELTAT * (A1(JU.K) + A2(JU,K)) 

TKM(JU) = THTEMP 

IF (K .EQ. KMID - IRATIO) TKM2(JU) = THTEMP 

GO TO 310 


GROUND SURFACE - INITIALLY 


fart | Lon eeOwaliGOe LOs2dO 
ZKM(J) = ZO(U,K) 


A4(JU,K) = -9.999D0 
AS(JU,K)..= -9.999D0 
A6(JU,K) = -9.999D0 


TH1(JU,K) = THO(JU,K) + DELTAT * COOLR 
TKM(JU) = THO(J,K) 


Ai(J,K) = -9.999D0 
A2(J,K) = -9.999D0 
GO TO 310 


GROUND SURFACE - SUBSEQUENT TIMES 


ZKM(JU) = ZO(J,K) 


A4(JU,K) = -9.999D0 
A5(U,K) = -9.999D0 
A6(J,K) = -9.999D0 


THTEMP = THO(J,K) 

THO(JU,K) = TH1(JU,K) 

TH1(U,K) = THTEMP + 2.D0 * DELTAT * COOLR 
TKM(JU) = THTEMP 


Ai(J.K) = -9.999D0 
A2(JU,K) = -9.999D0 
GO TO 310 


TOP BOUNDARY OF REGION 


A1i(JUJ,K) = 0.DO 
A2(U,K) = 0.DO 
A4(J,K) = 0.DO 
A5(uU,K) = 0.DO 
A6(JU,K) = 0.DO 
21 Gl Ke 40.00 


TH1(JU.K) = THO(UJ,K) 
GO pL0u3 40 
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ABOVE TROUGH OR RIDGE - INITIALLY 
230 IF? CITC) VEG. 4) 960 TO 260 
CALL ADVEC(ZO(JU;K + INC), ZO(U,K), ZKM(JU), DZ, WO(JU,K), VADV) 
A4(JU,K) = -VADV 
DIFZ = (Z0(JU,K - INC) - 2.D0*ZO(U.K) + ZO(JU,K + INC)) / DZ / 
1 DZ * KVZ 
DIFY = -2.D0 * ZO(U,K) * KHZ / DY / DY 
A5(JU,K) = DIFY + DIFZ 
A6(JUJ,K) = 0.DO 
Z1(JU,K) = ZO(JU,K) + DELTAT * (A4(JU,K) + A5(U,K)) 
ZKM(JU) = ZO(U,K) 
IF (K .EQ. KMID - IRATIO) ZKM2(JU) = ZO(U,K) 
CALL ADVEC(THO(JU,K + INC), THO(JU,K), TKM(JU), DZ, WO(U,K), 
1 VADV) 
DIFZ = (THO(U.K - INC) - 2.DO*THO(U,K) + THO(U,K + INC)) / DZ 
1 / DZ * KVT 
IF (J .GT. 1) GO’ TO 240 
DIFY = 2.D0 * (THO(U + 1,K) - THO(J,K)) / DY / DY * KHT 
GO TO 250 
240 DIFY = 2.DO * (THO(U - 1,K) - THO(JU,K)) / DY / DY * KHT 
250 Ai(U,K) = -VADV 
A2(JU,K) = DIFY + DIFZ 
THi(JU,K) = THO(JU,K) + DELTAT * (A1(U.K) + A2(U,K)) 
TKM(JU) = THO(J,K) 
IF (K .EQ. KMID - IRATIO) TKM2(JU) = THO(U,K) 
GO TO 310 
ABOVE TROUGH OR RIDGE - SUBSEQUENT TIMES 
260 CALL ADVEC(ZO(JU,K + INC), ZO(JU,K), ZKM(JU), DZ, WO(JU.K), VADV) 
A4(JU,K) = -VADV 
CALL DIFFOR(ZO(U,K + INC), ZO(JU,K), ZKM(JU), DZ, KVZ, DIFZ) 
TFEQUSE GT..PH) (GO: TO*270 
CALL DIFFOR(ZO(JU + 1,K), ZO(JU,K), -ZO(U + 1,K), DY, KHZ, DIFY) 
GO TO 280 
270 CALL DIFFOR(-ZTEMP, ZO(JU,K), ZTEMP, DY, KHZ, DIFY) 
280 ZTEMP = ZO(U,K) 
A5(U,K) = DIFY + DIFZ 
A6(JU,K) = 0.DO 
ZO(U,K) = 21(Jd,K) 
Z1(U,K) = ZTEMP + 2.DO * DELTAT * (A4(JU.K) + A5(U,K)) 
ZKM(J) = ZTEMP 
IF (K .EQ. KMID - IRATIO) ZKM2(JU) = ZTEMP 
CALL ADVEC(THO(U,K + INC), THO(JU,K), TKM(JU), DZ, WO(U,K), 
1 VADV) 
A1i(JU,K) = -VADV 
CALL DIFFOR(THO(U,K + INC), THO(JU,K), TKM(JU), DZ, KVT, DIFZ) 
TF? (0 £67.09) AG0 for29s0 
CALL DIFFOR(THO(JU + 1,K), THO(U,K), THO(U + 1,K}, DY, KHT, 
{ DIFY) 
GO TO 300 
290 CALL DIFFOR(THTEMP, THO(JU.K), THTEMP, DY, KHT, DIFY) 
300 A2(JU,K) = DIFY + DIFZ 


THTEMP = THO(UJ,K) 
THO(JU,K) = TH1(U,K) 
TH1(JU,K) = THTEMP + 2.DO * DELTAT * (A1(JU,K) + A2(J,K)) 
TKM(JU) = THTEMP 
IF (K .EQ. KMID - IRATIO) TKM2(JU) = THTEMP 
310 CONTINUE 
320 CONTINUE 
LFIEM ETTOTATPR. 2 ITTOT/ I TPR: MNE.*« OF0)3G0 TO 330 
IUNIT = 6 
WRITE (IUNIT,800) 
WRITE (IUNIT,700) TTOT 
CALL OUTPUT(IUNIT, A1) 
WRITE (IUNIT,710) TTOT 
CALL OUTPUT(IUNIT, A2) 
WRITE (IUNIT,690) TTOT 
CALL OUTPUT(IUNIT, TH1) 


Scr 


i 
(VOAY , dh, L708 .&C ‘ unmoee ke “be At. 
»e : 
\ 50 \ ((0Mn + A.LJos 4 OF ue 190. > (out Mah 


.(*,L)0W .fa .CL)MNT , 08 LIGHT LONE + F. he) 


wy vd (tut * UL OS~ , fe. L)os ,ta.) + LOS 807910 "a3 


TvN .S0 .(D)MNT O4.U)ONT , (owl + 8 )OHT)ROTITG JAD 


* 


Nite. Bb, 


tie. a) } 


LM. vor é 


08 
) sa Ago 


j 
7 OC i. € i 
= 


(C4, uaa © O84 pea) © rat sa +a 


ya.\ 40 \ > wet 


(Ww. LYO8 = (% )omne Martin - ab 


outs a LadeT + ¢a,uoNTsoGLe = Ci - 4b ON a} 
me We ‘ 
oss UY OB (tf . Te, : 
Hu * ¥O \ VO \ CCN. PONT = 08,9 4. 6)0HT? Sal wai90 
. $0 
¥ * yO \ ¥O)\. COMD)OMT ~ O80 = DIOHT) © Bele 
voayve *% 
S351G + ¥310 * 
(mM. L)CA * (807A) TATIIQ + (M.UJOMT = CH 
(4.G)OHT © 
(#, UJOHT = (LISMHT COLTART <=. QUA, O27. 
ore | 
r tupuoszeu2 ~ 3n0ia 20 neuont ‘~ 
“Aas 
LOW . £0 .CLIMNS | UNL 2OS . TOME + 8G 20RISSVGR Dae 


VOAV- © (UIRa 
Sv (S60 .(U}aNe , 08,0308  €2"T + 8 UPOS)ROIITE Jeeo 7 oe 
ors of o© (t . Te, LIM fi 


oss oT. - 

(yYaro .SH™ .YO ,SMyTS .(N.LIOS .SMaTS-)ROVIS ia 7 
(M,L)IOS @ SITS 

SRI + ¥R10 = (MLU)8 a . 

og.0 * (t,U)aa 

(w.UItS = (ALU - 

(¢%.U)24 © (8.U)8A) © TATIZQ * 00.8 © SMITS = (hu) 

aM27sS = (b ; 

site = (Gicwey (GTTAM? ~ GIMN .03. HP 

ca .(u)MNT , 08, )OHT . (OM) © &%,G20ONTIDaVOA 3380 

(VOAYV rm 

yoav- « (W.U)TA 

oes OF of (ft . Ta. MESSE) P 
wot » G2OMT . Ou JoMT Ja.t + G)onTyaotare a 

' ) a ae » 

oe OT 


VO .@GTHT . OM WPOMT SaNa TT IFO NTE rane of 
S910 + yaig = (A.u -* Wy 


(*,. VJORT = OM3THT > 
(x. L)THT © at 
¥Yicea * (H.UPTA) * TATIIC © O8.F + MITHT «= ren, — 
‘water © (yy 
MSTHT = (LUSMNT (TOTTFART - OM. . 02, 


wl 
“ 


o6t OT 00 (0.0 3M, AVTE\TOTTT + wake 1) 3b 


‘ (0O8, vse) atte 
‘yort (oor, frwut 7 a 
, “(ra . Ten) De B 4589 : 
> Oka , TIMI D, J4AS 
Torr toe J . 
ae riwMur d4AQ 


C 

Cc 

Cc 
330 


360 
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380 
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IUNIT = 8 

WRITE (IUNIT, 800) 
WRITE (IUNIT,770) TTOT 
CALL OUTPUT(IUNIT, A4) 
WRITE (IUNIT,780) TTOT 
CALL OUTPUT(IUNIT, A5) 
WRITE (IUNIT, 800) 
WRITE (IUNIT,720) TTOT 
CALL OUTPUT(IUNIT, A6) 


RELAX TO OBTAIN STREAMFUNCTIONS 


DO 380 I = 1, MAXITR 
ERROR = 0.DO 
DO 370 K = 2, KENDMI 
IF (K .GT. KMID .AND. (K - 1)/RATIO ~- (K - 1)/IRATIO .NE. 070) 
GO TO 370 
JFINAL = K + JCOR - 2 
IF (JFINAL .GT. JEND - 1) JFINAL = JEND - 1 
DO 360 J = 2, JFINAL 
IF (K .GE. KMID) GO TO 340 


DZ = DZ1 
INC ="1 
GO TO 350 
DZ = D2Z2 


INC = IRATIO 
PHIOLD = PHI(uU,K) 
RES = Z1(JU,K) + (PHI(U + 1,K) - 2.DO*PHI(U,K) + PHI(U - 1,K) 
)aA O¥e/- Dx 
RES = RES + (PHI(JU,K + INC) ~- 2.DO*PHI(U,K) + PHI(U,K - INC) 
)tfadzi 4 10z 
PHI(U,K) = PHIOLD + RES * RALPHA / (2.DO/DY/DY + 2.D0/DZ/DZ) 
PH = DABS(PHIOLD - PHI(uU,K)) 
IF (PH .LT. ERROR) GO TO 360 
ERROR = PH 
UM = J 
KM = K 
CONTINUE 
CONTINUE 
ITER = I 
IF (ERROR .LT. EPS) GO TO 390 
CONTINUE 
IF (ITTOT/TPR - ITTOT/ITPR .NE. 0.0) GO TO 400 
IUNIT = 9 
IF (IP1 .EQ. O) WRITE (IUNIT,.810) TTOT. ITER, ERROR 
IF (IP1 .EQ. 1) WRITE (IUNIT,800) 
IF (IP1 .EQ. 1) WRITE (IUNIT,740) TTOT, ITER, ERROR 
CALL OUTPUT(IUNIT, PHI) 


VELOCITY FIELDS ARE DERIVED FROM STREAMLINES 


DO 410 K = 1, KEND 
00° 440 Ut= 1) JEND 
VO(JU,K) = Vi(J,K) 
WO(JU,K) = W1(J,K) 
CONTINUE 
VM = 0.DO 
WM = 0.DO 


DO 510 K = 1, KEND 

IF (K .GT. KMID .AND. (K - 1)/RATIO - (K - 1)/IRATIO .NE. 0.0) 
, GO TO 510 

JFINAL = K + JCOR -' 1 

IF (UJFINAL .GE. JEND) JUFINAL = JEND 

DO 500 J = 1, JFINAL 
IF (K .EQ. KEND) GO TO 460 
TFAICK .EQ. 1) GO TO 470 
IF (K .EQ. J - JCOR + 1) GO TO 480 
IF (J .EQ. 1 .OR. JU .EQ. JEND) GO TO 440 
IF (K .GE. KMID) GO TO 420 
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DZ = D021 
INC = 1 
GO TO 430 
DZ = D022 
INC = IRATIO 
Vi(JU,K) = (PHI(JU,K + INC) - PHI(JU,K - INC)) / 2.00 / DZ 
W1i(U,K) = -(PHI(U + 1,K) - PHI(JU - 1,K)) / 2.DO / DY 
GO TO 490 
Vi(U,K) = 0.DO 
IF (J .EQ. JEND) GO TO 450 
Wi(JU,K) = -PHI(JUJ + 1,K) / DY 
GO TO 4980 
W1i(J,K) = PHI(JU - 1,K) / DY 
GO TO 480 
Vi(J,K) = 0.DO 
Wi(JU,K) = 0.DO 
GO TO 500 
Vi(JU,K) = (2.DO0*PHI(U,K + 1) - O.5DO*PHI(U,K + 2)) / 021 
IFS (JU .EQ. 1) Vi(dU,K) = 0.00 
Wi(JU,K) = 0.DO 
GO TO 490 
DPHIDY = (2.DO*PHI(U - 1,K) - O.5SDO*PHI(U - 2,K)) / DY 
DPHIDZ = (2.DO*PHI(JU,K + 1) - O.5DO*PHI(JU,K + 2)) / D021 
MAGNV = DSQRT(DPHIDY*DPHIDY + DPHIDZ*DPHIDZ) 
Vi(JU,K) = -MAGNV * DCOS(ALPHA) 
IF (J .EQ. JEND) V1(U,K) = 0.DO 
Wi(JU,K) = -MAGNV * DSIN( ALPHA) 
IF (J .EQ. JEND) W1(JU,K) = PHI(U - 1,K) / DY 
GO TO 4980 
IF (DABS(V1i(JU,K)) .GT. VM) VM = DABS(V1i(J,K)) 
IF (DABS(W1i(JU,K)) .GT. WM) WM = DABS(W1i(U,K)) 
CONTINUE 

CONTINUE 

IF (ITTOT/TPR - ITTOT/ITPR .NE. 0.0) GO TO 520 

IUNIT = 10 


IF (IP1 .EQ. O) WRITE (IUNIT,820) TTOT 
IF (IP1 .EQ. 1) WRITE (IUNIT,800) 

IF (IP1 .EQ. 1) WRITE (IUNIT,750) TTOT 
CALL OUTPUT(IUNIT, V1) 

TUNTT =" 14 

IF (IP1 .EQ. O) WRITE (IUNIT,830) TTOT 
IF (IP1 .EQ. 1) WRITE (IUNIT,800) 

IF (IP1 .E£Q. 1) WRITE (IUNIT,760) TTOT 
CALL OUTPUT(IUNIT, W1) 


VORTICITY AT SURFACE 


Z1C1% 2)" ="02D0 
Z1(UEND,KVLY) = 0.DO 
DO 560 K = 1, KVLYM1 
JBEGIN = 2 
IF (K .GE. 2) JBEGIN = K - 1 + JCOR 
JFINAL = K - 1 + JCOR 
DO 550 J = JBEGIN, JFINAL 
IF (J .LT. JFINAL) GO TO S30 
DWDY = (-1.5D0*W1(JU,K) + 2.DO0*W1(U - 1,K) - O.5D0*W1i(uU - 2,K)) 


DwOY = 0.DO 
zi =4 7 5D0#VICU,.K) + 2.00*V1I(U Ko +~1) -~O. 5D0*V1(dU.K +°2)) 


Z1C ORK) 
Z1(JU,K) DWDY - DVDZ 

CONTINUE 
CONTINUE 
IF (ITTOT/TPR - ITTOT/ITPR .NE. 0.0) GO TO 570 
IUNIT = 8 
WRITE (IUNIT,730) TTOT 
CALL OUTPUT(IUNIT, Z1) 
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640 
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780 
790 
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810 


820 
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850 
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EEeGLP } bO~ 


gf Po as | 
GO TO 590 
IP? =.0 


ERION SATISFIED ? 


-1) GO TO 580 


TE. (VM, .EQ.. .0.0 .OR..WM .EQ. 0.0) GO TO 640 


IF (VM .GT. 
DEL UY sabes 
OEE TZS=.D 24 
DTOLD = DELT 
LET COELTLALs. 


DELTAT = DEL 


4.0 .OR. WM .GT. 4.0) GO TO 840 
VM / 2.D0 * 0.800 

/ WM / 2.D0 * 0.8D0 

AT 


iis DELLY o AND. DECTAR ULI. DEETZ) GO) TO 620 


Ales 2a 


WRITE (6,610) DELTAT, TTOT 


FORMAT (’ “, 

IF (DELTAT . 

GO TO 600 

IF (DELTAT . 

CiCoasO 

DO 630 K = 1 

DO 630 J = 

THO(UJU,K) 
ZO0(JU,K) 
vo(J,K) 
WO(UJU,K) 

CONT INUE 

Le CT TOT, “GE 

GO TO 160 

i Wage | 

IF (TTOT .GE 

GO TO 160 

FORMAT (/’ “, 


1EG K)’) 


FORMAT 
FORMAT 


~N 
x< 


~~ O- 


FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 


are are) - 3 


S hes Bos Gee S85 
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FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 
FORMAT 


FORMAT 
FORMAT 
IUNIT = 7 
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‘DELTA T REDUCED TO’, F6.2, ‘AT TIME =’, F10.2) 
LT. 1.D0) GO TO 840 
EQ. DTOLD) GO TO 640 
, KEND 
1, JEND 
= TH1(J,K) 
7h) C) EO) 
= Vi(J,K) 
= W1(J,K) 
. TEND) GO TO 840 
. TEND) GO TO 840 
‘HEIGHT(M. ) PRES(MB. ) TEMP(DEG kK) POT TEMP(D 
ES 4.40% F510 (GX A171) 
“AFTER”, IS, “ITER., MAX. ERROR =’. “F12-10. 
OCCURS Ale cK Sol Stag ee Lae eee) 
‘POTENTIAL TEMPERATURE (THETA) AT TIME = ’, F7.2) 
“ADVECTION TERMS (THETA) AT TIME = ‘’, F7.2) 
‘DIFFUSION TERMS (THETA) AT TIME = ’, F7.2) 
"BUOYANCY TERM (Z) AT TIME = ’, F7.2) 
"'YORTICITY, (2) AT TIME: =~, F722) 
‘STREAMFUNCTION (PHI) AT TIME = ’, F7.2, ’ (AFTER ’ 
ITERATIONS. ERROR =. fF 12-10, “)7) 
“CROSS-VALLEY WIND SPEED (V) AT TIME = ‘’, F7.2) 
“VERTICAL WIND SPEED (W) AT TIME = ’, F7.2) 
‘ADVECTION TERMS (Z) AT TIME = ’, F7.2) 
‘DIFFUSION TERMS (Z) AT TIME = ’, F7.2) 
) 
’STREAMFUNCTION (PHI) AT TIME = ’, F7.2, ’ (AFTER ’, 
ITERATIONS, ERROR =’, F12.10, ’)’) 
“CROSS-VALLEY WIND SPEED (V) AT TIME = ’, F7.2) 
“VERTICAL WIND SPEED (W) AT TIME = ’, F7.2) 


WRITE (IUNIT,850) ITC, TTOT, DELTAT 


FORMAT (I1, 
IF (ITC .EQ. 


2.7.2) 
1) GO TO 860 


WRITE (IUNIT,690) TTOT 

CALL SAVE(IUNIT, THO) 

WRITE (IUNIT,730) TTOT 

CALL SAVE(IUNIT, ZO) 

WRITE (IUNIT,810) TTOT, ITER, ERROR 

CALL SAVE(IUNIT, PHI) 

WRITE (IUNIT,820) TTOT ° 
CALL .SAVE(IUNIT, .VO) 

WRITE (IUNIT,830) TTOT 

CALL SAVE(IUNIT, WO) 


GO TO 870 


maT TOO 8§=6 (4 -WS0PGMET ( aMj23n4 {. 4) THOTRM’ a) 


oss of da (T3390 . TH. rarsag oe 


ToTT Th tl 
(¢.o1r9 ,** SMIY TA’ .$.89 ,°O? Gdaouoga - AT, 
O88 OF OD (0G,) . 


ob® OT Oa \as0T@ .08 airy 
wt 4 


uae 2 00 
OnaL .t. i 
(A, LIPHE # 
O, &)tS ar 
(a, une a Bs ay 
(MWe 2. is 
~~ Th 


One oT ce kOwsT . 30, 


oh8 GOT DD (aver 38. 


CAN. 083, XOIS .hvG «KOR 5 PMA xs mpl OF 
OM .St3 ,° « SOSA KOM , ARTE’ (C1 ,'RBPRAT AF Yanga se 


(*¢* ,@7 .*.' 96h 2%) 2 ta, TA 29990 OA ~ 


nye deal 
(s.t4 .° | SRIT TA (ATI4T) SMIFASSaeeeT JALTMRTOR ony : 
(¢.T3 ,° = SMIT TA (ATSMT) J2WRST WOITOSVGA’ .*.> 
(s.t3 ,* = 3MIT TA (ATSMT) 2MR3T VOTZURIIC™ th 


(¢.¥2 ,.* = SMES TA (5) XILOPTROV' 
RITIA) * .$.09 .* = SMITE TA (IMS) MOI TOMUAMASATS* 


° 
’ 
(¢.t9 ,* = SMIT Ta (5) MAST YOMAYOUB? ,* 
Ld 
® 
('¢* ,OF. era ,* © RORMS.2MOTTASaTI * 


(c.3 <= 3MiT TA (VV) GF 1392 OMI VY3JJaV-22089' , 
(¢.T3 .* @ SMIT TA OW) 33 392 CYIw Js arteasv' . 
~ $8.99 2 WIT TA (S} @MA3ST MOITOSVGA® . 
(¢.T3 .° = BMIT TA (CS) BmHST WOT2UNAIO’., 
Q3T%A) * .¢.79 |" © BMIT TA CLN@) KOTTOMUIMAaRTe* tp 
(‘(* (OF Shay, “= RORRD, SMODTARSTA © .E% ‘te 
(¢ 2 .° © SMYT TA (¥) GaaRe OMIW VS JjAN-220R0' . ‘7 * ie is 
(o'83 .) = aMfT TA (a) S892 GnW JADEtRaV’ P+) ale 
v= Th 


JATIIN . TOTY .aTt (088, Tiwi) at 
i¢.t3c) ta) 4 
@oe of DD (? .05. 
jot? (os7,TiMutD ” 
(OHT ees JAD 

TOTT (Ott. Tru T1e “ie 

(OS, , Trwud javae. 

noers ,Aarh. OTT (ors, Taut). 4 i 

(THA TAMUE WAR. Naa 
‘Fort (Oss, TIMUE) # re 7 
{ov . TIAUE VAR i 7 
“tor (oEe, Yiwu) 5 
tow wile) 7 


| ie 7 7 


aan 


aaa 


omen) 


860 


870 


10 


10 


20 


30 
40 
50 


60 
70 


™ = TTOT - DTOLD ~ 
WRITE (IUNIT,690) TM 
CALL SAVE(IUNIT, THO) 
WRITE (IUNIT,690) TTOT 
CALL SAVE(IUNIT, TH1) 
WRITE (IUNIT,730) TM 
CALL SAVE(IUNIT, ZO) 
WRITE (IUNIT,730) TTOT 
CALL SAVE(IUNIT, 21) 


WRITE (IUNIT,810) TTOT, 


CALL SAVE(IUNIT, PHI) 
WRITE (IUNIT,820) TM 
CALL SAVE(IUNIT, VO) 
WRITE (IUNIT,820) TTOT 
CALL SAVE(IUNIT, V1) 
WRITE (IUNIT,830) TM 
CALL SAVE(IUNIT, WO) 
WRITE (IUNIT,830) TTOT 
CALL SAVE(IUNIT, W1) 
STOP 

END 

SUBROUTINE ADVEC(ABO, 


FORWARD-TIME, UPSTREAM-SPACE FORMULATION OF ADVECTION TERMS 


AMID, 


ITER, 


TMPLICIT REAL *S(A’ - H.O- 2) 
lGeGvVels.Gks ,0:00): GO.L0, 10 
RES = VEL * (ABO - AMID) / D 


RETURN 


RES = VEL * (AMID - ABEL) / D 


RETURN ’ 
END 
SUBROUTINE DIFFOR(ABO, 


FORWARD-TIME, CENTERED-SPACE FORMULATION OF DIFFUSION TERMS 


AMID, 


IMPLICIT REAL*8(A - H,O - Z) 


R= AK /ODO/OD 


ABEL, 


ERROR 


ABEL, 


RES = R * (ABEL - 2.DO*AMID + ABO) 


RETURN 
END 


SUBROUTINE OUTPUT(IUNIT, 


DATA IS WRITTEN ONTO VARIOUS LOGICAL UNITS 


A) 


IMPLICIT REAL*8(A - H,O - 2) 


DIMENSION A(16,28), FOR1(10), 
COMMON TTOT, FOR1, FOR2, 


KK = KEND + 1 - KL 
Jie =— VIB 

Jie —eON OI +d 
J3e=suU4 +edUJB 

J2 = J3 


IF (J2 .GT. JEND) U2 = 


KL, 


JEND 


FOR2( 10) 
JEND, UCOR, 


IF (Ji .NE>S 1) WRITE (CIUNIT.20) 


FORMAT (/17) 
DO 70 K = KK, KEND 

= KEND + 1- K 
IF (U1 .GT. (KBACK - 
ire CU2eGie (KBACKS— 


Tre tIGNIT, EQ. 6).-GO" T0'-30 


WRITE (IUNIT,FOR1) (A(U,KBACK) ,J=JU1,J2) 


GO TO 50 


1 + JCOR)) GO TO 70 
1 + JCOR)) u2 


TF CAG) Ghee OO) GU ETLOL4O 


WRITE (IUNIT,FOR1) (A(JU,KBACK) , J=JU1,UJU2) 


GO TO 50 


WRITE (IUNIT,FOR2) (A(U,KBACK) , J=JU1,U2) 
IF (KBACK .GT. KMID) K 
IF (KBACK .GT. KMID) WRITE (IUNIT,60) 


FORMAT (’ ’, /) 
CONT INUE 


VLG A 1 


LE (Jae JLT... JENDpe GOTO S10 


RETURN 
END 


D ’ 
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SUBROUTINE SAVE(IUNIT, A) 


DATA IS WRITTEN ONTO LOGICAL UNIT 7 FOR CONTINUATION OF RUN 
AT A LATER TIME 


IMPLICIT REAL*8(A - H,O - 2Z) 
DIMENSION A( 16,28), FOR1(10), FOR2(10) 
COMMON TTOT, FOR1, FOR2, KL, JEND, JCOR, KEND, KMID, UJB 
KK = KEND + 1 - KL 
set 
DO 20 K KK, KEND 
KBACK KEND + 1 - K 
J2 = KBACK - 1 + JCOR 
IF (J2 .GT. JEND) J2 = JEND 
WRITE (IUNIT,10) (A(JU,KBACK) , J=JU1,J2) 
FORMAT (1X, 16D23.15) 
IF (KBACK .GT. KMID) K = K + 2 
CONTINUE 
RETURN 
END 
SUBROUTINE UNSAVE(IUNIT, A) 


DATA IS READ FROM LOGICAL UNIT 12 FOR CONTINUATION OF RUN 
AT A LATER TIME 


IMPLICIT REAL*8(A - H,O - Z) 
DIMENSION A(16,28), FOR1(10), FOR2( 10) 
COMMON TTOT, FOR1, FOR2, KL, JEND, JCOR, KEND, KMID, JuB 
KK = KEND + 1 - KL . 
Ji = 1 
DO 20 K KK, KEND 
KBACK KEND + 1 - K 
J2 = KBACK - 1 + JCOR 
IF (J2 .GT. JEND) J2 = JEND 
READ (IUNIT,10) (A(JU,KBACK) ,J=U1.J2) 
FORMAT (1X, 16023. 15) 
IF (KBACK .GT. KMID) K = K + 2 
CONTINUE 
RETURN 
END 
SUBROUTINE SUBK(JEND, KEND, UCOR, KVLY) 


DEFINES DIFFUSIVITIES (ALTER AS DESIRED) 
IF K = K(Z), HOWEVER, GOVERNING EQUATIONS MAY CHANGE 


IMPLICIT REAL*8(A - H,O - Z) 
REAL*8 KHZ, KVZ, KHT, KVT 
COMMON /A/ KHZ, KVZ, KHT, KVT 


KHZ = 1.D0 
KA =s. DO 
KVZ = 0.DO 
KVT = 0.DO 
RETURN 


END 
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